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METHOD OF FAN SOUND MODE STRUCTURE DETERMINATION 
COMPUTER PROGRAM USER'S MANUAL 
MODAL CALCULATION PROGRAM 


by 

G. F. Pickett, R. A. Wells end R. A. Love 
Pratt & Whitney Aircraft Group 

1.0 SUMMARY 

This computer user’s manual describes the operation and the essential features of the Modal 
Calculation Program, the second of two programs developed under the Method of Fan Sound 
Structure Determination Program, NAS3-20047. Jointly the two programs are used to deter- 
mine the coherent modal structures of inlet sound fields. The purpose of the Modal Calcula- 
tion Program i* to calculate the amplitude and phase of modal structures by means of acoustic 
pressure measurements obtained from microphoes placed at selected locations within the fan 
inlet duct. These locations are determined by the first of the two programs. In addition, the 
Modal Calculation Program also calculates the first-order errors in the modal coefficients that 
are due to tolerances in microphone location coordinates and inaccuracies in the acoustic pres- 
sure measurements. 


2.0 INTRODUCTION 

New fan designs for modem high bypass ratio commercial engines utilize blade-vane interac- 
tion theory to the extent possible for contolling the propaga ion of interaction noise. Cur- 
rently. this theory defines the modes that can propagate, but has not been developed to the 
extent that it can reliably predict the strengths of the propagating modes. 

Further noise reduction could be achieved if the propagating modal structure were quantified. 
Once the modal structure were defined, an analytical system for acoustic-treatment design 
could be utilized to optimize treatment for a given modal structure, to produce more efficient 
schemes. In addition, the modal structure could be employed to verify developing theories of 
fan noise generation. To provide this capability by means of measured data the Method of 
Fan Sound Mode Structure Determination Program (NAS3-20047) was undertaken. The 
method would be utilized until a valid fan noise generation model on a model basis becomes 
available. 

The theory upon which fan spinning mode theory is founded was presented in 1961 by Tyler 
and Sofrin (ref. I ), following extensive analytical and experimental studies. Later, Sofrin and 
McCann (ref. 2) derived the general form of a coherent acoustic wave in an infinitely long cy- 
lindrical duct which extended the theory to include effects of axial flow. This equation 
expresses the coherent acoustic pressure at locations in the duct as a function of the ampli- 
tude and phase of the propagating modes comprising the sound field. These purely coherent 
signals, which are due to the contributions of the constituent modes, are extracted from the 
overall signal by enhancement techniques adapted at Pratt & Whitney Aircraft - the advan- 
tages of utilizing signal enhancement is discussed by Posey in reference 3. 


Both the analytical expression derived for a general coherent acoustic wave and a signal en- 
hancement technique form the basis for developing a method tp determine fan sound mode 
structures. The method, in principle, is capable of determining the amplitude and phase of 
all modes that can propagate at a given frequency. In practice, the number of modes that 
can be determined is limited by the storage capacity and the running time of the computer 
and by measurement and location accuracy. 


The method for determining fan sound mode structure (ref. 4) requires two computer pro- 
grams: a Microphone Location Program (MLP) and a Modal Calculation Program (MCP). This 
User’s Manual describes the MCP; the MLP is presented in a companion Manual. 

The MLP identifies microphone locations in the duct for measuring acoustic pressures for 
input to the MCP that will insure a numerically stable solution. The MCP calculates modal 
structures from acoustic pressure measurements and calculates coefficients that can be used 
to determine the sensitivity of the modal calculation procedure to first-order errors in 
acoustic pressure measurements and microphone placement. 

In the following sections, the algorithm for the modal calculations and the program elements - 
such as subroutines, functional elements, and principal element interrelationships - are di* 
cussed. A description of the input parameters is included. The output format is also described 
and illustrated by a sample case. Finally, a listing of the program code is provided in 
Appendix B. 
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3.0 PROGRAM DESCRIPTION 


3.1 ALGORITHM 


The Modal Calculation Program is an algorithm for calculating the modal structure from in- 
put data composing acoustic pressure measurements and a finite set of modes. The general 
form of any coherent acoustic wave in an infinitely long cylinderical duct having uniform 
axial flow can be written as the real part of 


P(x,r,0:t) = SA m ^Efk^r) e 
1 iniie Set 
of Modes 



and 

Kx= Mx (“ /c > 1 \/ (^c) 2 ' ( 1 - Mx 2 ) k ^ 2 
‘-Mx 2 

where the notation is consistent with reference 1 . 

Equation 1 can be written in matrix form where the measured pressures are obtained from 
microphone locations identified by the MLP. The equation system is solved in the usual 
manner by matrix inversion. The output from this procedure is the amplitude and phase of 
the coherent acoustic duct modes comprising the inlet sound field. 

The input to the program consists of: the sound field in the duct comprising N acoustic duct 
modes, the geometric parameters (e. g. duct radius, hub-tip ratio), test parameters (e. g. freq- 
uency, axial Mach number, speed of sound), and measured acoustic pressure amplitude and 
phase at locations identified by the MLP. The characteristic numbers that include the eigen 
value k' p , the axial wave number kx, and the value of the eigen function E (k^ r) 
are calculated. 

In addition, this equation requires the input of acoustic pressure measurements, the number 
cf which exactly equal the number of specific modes. A set of equations can then be estab- 
lished with the number of equations equaling the number of acoustic measurements. This set 
was written in matrix form with the matrix coefficients a function of the particular modes 
comprising the sound field and the microphone locations. 

If the determinant of the equation system is non-zero, a set of independent equations exists. 
This equation system in principle can be inverted in the usual way to solve for the unknown 
amplitude and phase of the particular modes comprising the sound field. A Gaussian elimina- 
tion procedure is used to reduce the equation system to a triangularized matrix for solution 
of the complex modal coefficients. The overall pressure at any location in the duct can be 
calculated from the information in the modal structure. 
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Once the modal structure has been determined, a set of influence coefficients (ref. 4) is cal- 
culated. These coefficients can be used to determine the errors in modal amplitudes and 
phases that are the results of first-order inaccuracies in measured pressures and the tolerances 
in microphone placement. 

As an option, the MCP can also calculate the resultant sound field at any specified duct loca- 
tion based on - 3 iven modal structure. This modal structure is supplied by the user either 
arbitrarily or as output from an analytical prediction deck. 

3.2 PROGRAM OVERVIEW 

The Modal Calculation Program comprises six major sections which are utilized in part or 
whole to accomplish the objectives of the two possible modes of operation. These six major 
sections are: 

1) Input - The input of all data is by the NAMELIST specification, and the internal para- 
meters are initiated for program execution. 

2) Characteristic Number Calculation - The characteristic numbers K m °^ and 
are calculated using the procedure described in Appendix A. 

3) Mode Amplitude and Phase Calculation - The coherent acoustic wave equation system,, 
e.g. (1), in an infinitely long cylinderical duct with uniform axial flow is solved using a 
Gaussian elimination procedure for the modal amplitude and phase. 

4) Sensitivity Coefficient Calculation - Standard deviations due to the first-order independ- 
ent errors in the measurement of both the acoustic pressures and the microphone coordin- 
ates are obtained for the error in the modal amplitudes and phases. 

5) Overall Pressure Calculation - Resultant pressure amplitude and phase are calculated at 
the desired prediction locations using the amplitude and phase of the constituent modes 
comprising the sound field. 

6) Output - All results from the program calculations are printed. 

The interrelationships between the six major sections and their utility for each option is il- 
lustrated in Figure 1 . As input, both options require a specific mode group, inlet geometry, 
and test condition to calculate characteristic numbers. One option, “A”, requires additional 
input in the form of acoustic pressure signals at selected duct locations to calculate the modal 
structure comprising the sound field. Additionally, influence coefficients, which are functions 
of the modal structure, are calculated. The other option, “B”, requires that the modal struc- 
ture be specified as input. In both options, the amplitude and phase of the constituent modes 
are utilized to calculate the overall acoustic pressure at any duct location. The results from 
both options are printed by the output section. 


4 



t-i! 


; > 
i 



'j 


r-f 


j 

4 




•i 


• J 

r.\ 

v 







I 


i 


i 


MODE AMPLITUDE 
AND PHASE 
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Figure 1 Program Overview 


3.3 PROGRAM SUBROUTINES AND FUNCTIONS DESCRIPTION 


The subroutines and functions used in the six program sections presented in Section 3.2 are 
listed belcw; the purpose of each subroutine or function is described. Also as appropriate, 
principal-element diagrams of the more complicated sections are presented and discussed. 


i 
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Input Section 

The NAMELIST format is used to input data for execution of the computer program. This 
form of input is described in Section 3.4.1. The input variable names are listed in Section 
3.4.2, including a description of their purpose. All input is read into the program by the 
following subroutine: 

INPUT - This subroutine inputs data for each case and sets up the necessary in- 

ternal parameters. 

Characteristic Number Calculation Section 


Expressions are derived in Appendix A for solving two simultaneous equations that define the 
characteristic numbeis k and Q » . A principal-element diagram is presented in Figure 2 
to illustrate the functional elements mat lead to a determination of these numbers. Initially, 
the order of the Bessel functions is determined from the circumferential order of a particular 
mode. The J m and Y m Bessel functions are evaluated, as appropriate, depending on the value 
of the duct hub-tip ratio. Finally, the characteristic numbers are calculated by solving the 
simultaneous equations comprising the Bessel functions. The subroutines and functions 
utilized in this section are: 

KQCAL - This subroutine calculates the characteristic numbers K^ and Q 


KMUCAL 

EMUCAL 


This subroutine is used by KQCAL to calculate the characteristic num- 
ber 1^ m /i • 


This subroutine calculates characteristic E-function values fora particu- 
lar radial value, r ' - r/b. 


FALZIP 

BESL1 


BESL2 


BESJ 

BESY 


This function solves for a root of a given function using a combination 
of false position and bisection techniques 

This function is used by KMUCAL to calculate values of K ^ for the 
equation which defines the system of differential equations. 

”tfr ^m ^ m/i^ + ^ mp^ r ' l Y m ^m/a^ = ° 

4 + <>w7r' 1Y "' (,K W 11 =0 

for a hub-tip ratio not equal to zero. 

This function is used by KMUCAL to calculate values of K^ for the 
equation which defines the above system of differential equations for 
a hub-tip ratio equal to zero. 

This subroutine calculates values of the Bessel function of the first kind. 
This subroutine calculates values of the Bessel function of the second kind. 
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Mode Amplitude and Phase Calculation 

The modal amplitude and phase are solved by matrix inversion techniques from data that in- 
cludes pressure measurements at selected microphone locations. The equations that define 
the matrix coefficients and a description of the procedure for fan sound mode determination 
was presented in Section 3.1 - Algorithm. To illustrate the functional elements that lead to 
a solution of the modal coefficients, a principal-element diagram is presented in Figure 3. 
Initially, the matrix coefficients, which are functions of the particular modes comprising 
the sound field and the microphone locations, are calculated. This equation system is solved 
by a Gaussian elimination method for the modal coefficients. The mode amplitude and 
phase are then extracted from these complex pressure vectors. The subroutines used in the 
calculation procedure are: 

S0LVF - This subroutine set ups and using SIMKCQ solves the acoustic wave equa- 
tion matrix for the modal amplitude and phase. 

S1MECQ This subroutine solves a N x N system of simultaneous equations having 

complex coefficients, using a Gaussian elimination method. 

Sensitivity Coefficient Calculation 

The Sensitivity Coefficient Calculation procedure is illustrated in the principal-element dia- 
gram presented in Figure 4. 

An important element in this procedure is the calculation of influence coefficients, which re- 
flect the sensitivity of mode amplitude and phase calculations to first-order errors in pressure 
measurements and microphone placement - the derivation of the influence coefficient is pro- 
vided in reference 4, Section 3.4. 

Because the inverse-matrix element is a common term in each expression, the procedure is 
initiated by calculating the inverse matrix. The influence coefficients are calculated next 
as a function of the modal structure and press) ire measurements. The specific error in the 
modal amplitude and phase due to one of the five possible measurement errors is calculated 
from the product of the error in the measured quality and the root-sum-square of the in- 
fluence coefficients. Finally, the error in a particular mode amplitude and phase is obtained 
as the combined effect of each measurement error. 

The subroutines utilized in the Sensitivity Coefficient Calculation procedure are: 

SENSTY - This subroutine calculates the standard deviations of the modal ampli- 

tude and phase for errors associated with pressure measurement and 
microphone location. 

INVERT - This subroutine inverts a complex N x N matrix. 
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Figure 3 Principal- F lenient Diagram Mode Amplitude and Phase Calculation 
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Figure 4 


Principal-Element Diagram - Sensitivity Coefficient Calculation 



Overall Pressure Calculation 


The overall pressure at any location in the duct is obtained from the modal structure. The 
procedure for overall pressure calculation is illustrated in the principal-element diagram shown 
in Figure S. The procedure summarizes the pressure contribution of each mode at a location 
in a duct defined by the user. The resultant amplitude and phase are then extracted from the 
complex pressure vector. Since this calculation is performed in the MAIN routine there are no 
subroutines or functions to list. 

Output Section 

The output format and the variables from the Modal Calculation Program are discussed in 
Section 3.5.1 and a sample case for three propagating modes is providing in Section 3.5.2. 

Both Sections 3.5. 1 and 3.5.2 address the two possible modes of operation that can be exe- 
cuted with the program. Results from the computations are printed by the subroutines 
listed below after all angles are converted to within the range of 0° to 360° . 

PRINT - This subroutine prints input and resultant values. 

ANGPQS - This subroutine converts negative angles to positive angles in the range 

0° to 360° for printing. 
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Figure 5 Principal-Element Diagram - Overall Pressure Calculation 
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3.4 INPUT DESCRIPTION 






3.4.1 Input Format 

The NAMELIST format is used to input data into the Modal Calculation Program and con- 
sists of a list of parameter names grouped under an identifying name: &INDATA. The 
parameter names correspond to variables - single variables and matrix elements - used in 
the program. These variables are set by specifying both the parameter name and its value. 

A feature of this type of input is that all associated parameters need not be specified. Any 
parameter not specified in the input retains its value from the preceeding case or the default 
value if the input is for the first case. 

NAMELIST input for each case is identified by the characters &INDATA in Columns 2-7 
of the first input card. Beginning in Column 9, parameters may be set using the format: 

Parameter Name = Constant 

The constant may be either a real or integer value and must be followed immediately by a 
comma. Parameter names, assigned values, or necessary commas must not extend beyond 
Column 72; and names of values cannot be continued on a subsequent card. Embedded 
blanks are not permitted in either the parameter name or constant value. Parameter names 
and their associated values may be specified in any order. The characters &END signify the 
end of the input for a particular case. If additional cards are required, parameters names 
must begin in Column 2. 

A sample of this form of input for three microphones is presented in Figure 6. 

3.4.2 Input Parameters 

A sign convention was adopted for assigning positive or negative values to the input para- 
meters. Any input parameter not addressed in this discussion is a positive value. The sign 
convention is formulated with respect to a cylindrical coordinate system that is consistent 
with the derivation of the coherent acoustic wave propagation model, its unit vectors are 
designated by the directions: axial - x. circumferential - 6, radial - r. 

A cons at radius, annular duct is aligned with respect to this coordinate system in such a 
way that the positive axial unit vector is in a direction opposite to the flow. Thus, the Mach 
number of a uniform axial flow is always designated by a negative value to denote the axial 
Uuw rate in the negative axial direction. A positive circumferential unit vector projects in 
the direction that the rotor spins, and a negative vector projects in the counterrotating direc- 
tion. Finally, the radial axis projects perpendicular to the centerline of the duct; thus, radial 
values are positive. 

Each mode is characterized by three parameters which represent the circumferential and rad- 
ial pressure distribution and its propagation direction. A specific mode is uniquely defined 
by the parenthetical notation (M, n). The M defines a periodic circumferential pressure dis- 
tribution with M number of lobes. Positive integers represent a corotating M-circumferential 

13 




7 


T r 

f 


* 9 


i 






lobe pattern with respect to the rotor direction, and negative M integers refer to counterrotat- 
ing modes. The radial mode index p corresponds to the radial pressure distribution. These 
values are always non-negative integer numbers with high integer values indicating large pres- 
sure variations with respect to the radius. 

The modal propagation direction in an inlet or discharge duct can be either an incident wave 
pn legating from the fan or a reflected wave propagating towards the fan. Wave propagation 
in , moving medium is similarly effected by the flow rate for modes that are propagating with 
. <1 against the flow direction. Hence, the input variable IDIR designates wave propagation 
v th respect to the flow direction. Positive values denote waves propagating in the opposite 
direction with respect to the flow such as incident waves in the inlet duct and reflected waves 
in the discharge duct. Modes that propagate in the same direction as the flow are designated 
by a negative value for the input parameter IDIR. 

A; signing of values to the input parameters will now be considered. 

Si ice a determinative equation system is required, the number of mode indices, wave direc- 
tion indicators, microphone coordinates, and measured pressures must be equal. When option 
B is utilized, the number of mode indices, wave direction indicators, and modal amplitude and 
phase values must correlate. These input parameters are listed in several tables at the end of 
th s section. Each parameter has a corresponding description that is sufficient for assigning 
a value to these input parameters. However, assigning a value to the coefficient parameters 
for the standard deviation in measurement errors is not as straight forward as the previous 
parameters. The following discussion is provided to assist the user when assigning values to 
these variables. 

The deviation coefficients for microphone location errors r:e the tolerances in the three co- 
ordinates: axial - x, radial - r, and circumferential - 6 . These errors are related to the tolerance 
of a measurement - such as a micromenter - for determining the location of a microphone. 
Specifically, a user can estimate the microphone location standard deviation by assuming a 
high confiden level - such as ninety-five percent - to be associated with the number of 
significant digits used to define the pressure measuring coordinates. The standard deviation 
coefficients can then be computed from this information. For example, if a 95 percent con- 
fidence level is assigned to an axial measurement accuracy of 0.005 centimeter, the standard 
deviation (68.3 percent confidence level) is about 2.5 x 10 - ^ centimeter. 

The error deviation coefficients for acoustic pressures include the two components amplitude 
and phase which correspond to the measured resultant pressure at any duct location. Two 
mechanisms can gmente errors that affect the measurement of resultant pressure. One type 
of error is d. . to both response characteristics of the measuring device and repeatability of 
the cohe^at signal. The second type of error is caused by measuring contributions from 
mode . not included in the calculation for determining the modal structure. A user can esti- 
ma„ the former pressure measurement error in a similar manner as previously presented for 
it.icrophone location measurement errors. 
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A standaid deviation can be computed by assuming a high confidence level to be associated 
with the combined inaccuracy of both pressure amplitude calibration errors and an t'rrcr 
attributed to the repeatability of enhanced pressure signals during a period of time. In prac- 
tice, however, this category of errors is small and can be minimized by requiring reasonable 
experimental procedures. 

The second mechan. m that can generate pressure measurement errors was not encountered 
in the previous category of location measurement errors. Ideally, the contribution from 
modes that are unlikely to control the duct sound field will not hinder the determination of 
fan sound mode structures. In practice, however, these modes have to be anticipated and their 
impact quantified if a meaningful standard deviation for the modal coefficients is to be cal- 
culated. This mechanism, which can be perceived as a measured pressure error, is difficult to 
assess prior to an experimental program. A general expression for this standard deviation is 
presented in Appendix E of reference 4. The actual value for the standard deviation used 
as input to the modal calculation program should be obtained from that general expression. 

A description of the input variables for operating the Modal Calculation Program is provided 
in Tables I, II, and III: Table I - General Parameters; Table II - Test Geometry and Condi- 
tion Parameters; Table Ill-Error Deviation Coefficient Parameters. Under the column head- 
ing “Variable Type”: th ' letter “R” indicates that the number is real and contains a decimal 
point; the letter “I” indicates the number is an integer and does not have a decimal point. 
“Default Values” are also delineated and indicate the value of the parameter that is internally 
initialized prior to the program execution. Parameters not specified in the input for the first 
case retain this value. Although the default values are expressed in units of the English System, 
the computer program can be executed with data in any consistent system of units. 

3.5 OUTPUT DESCRIPTION 


3.5. 1 Output Format 

The output from the Modal Calculation Program is organized into four sections: Input 
Variables, Modal Amplitude and Phase Calculation, Sensitivity Coefficient Calculation, and 
Characteristic E-Punction Values. All four sections are included as output when either option 
is requested by the input. The printout for a sample case is provided in Appendix C to 
illustrate the output format. 

The Input Variable Section includes the value of the various parameters supplied by the user. 
The parameters that define the modal structure - the circumferential and radial order, and the 
wave-direction indicator - are listed. The reference pressure for converting the modal amp- 
litude and resultant amplitude to decibels is also output in this section. The test geometry 
and conditions subsection lists various parameters that define the fan duct geometry and 
operating conditions observed during the experimental program. These parameters include 
the duct radius, duct hub-tip ratio, axial Mach number, and frequency. 

The Modal Amplitude and Phase Calculation section includes both parameters that were pro- 
vided by the user and the results from the calculation procedure. In this section, the user 
obtains the modal amplitude in units of pressure and decibels and the modal phase in units 


ORIGINAL PAGE B 
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TABLE I 


GENERAL INPUT PARAMETERS 


Input 

Name 


LOCM 


LOOP 


1EMU 


PREF 


X0 a) 

THO a) 


M(l) 

M(2) 

M(3) 


Variable 

Oefault 


Type 

Value 

Description 

l 

2 

Number of microphones or modes. (Less than 
or equal to fifty). 

1 


Number of prediction locations. (Less than 
or equal to fifty). 

1 

0 

Print indicator for characteristic E-function 
value. 



0 * No print 

1 = Prim 

R 

2.9 x 10' 9 

Reference pressure to convert pressure to 
decibels. 

R 

0.0 

Axial coordinate of the reference location. 

R 

0.0 

Circumferential coordinate of the reference 
location, (degrees) 

l 

-2 

Circumferential mode index. 


0 

(Input NLOC values) 


M(50) 

MUS(I) I 

MUS(2) 

MUS(3) 


0 

0 Radial mode index. 

1 (Input NLOC values) 
0 


MUS(50) 

ID1R(1) 1 

IDIR(2) 

IDIR(3) 


1D1R(50) 


Note a): Final character is a *eio. 


0 

I Mode propagation direction indicator. 

1 (Input NLOC values) 

0 

1 = opposite flow direction 
1 = with flow direction 

0 
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TABLE II 


* 


TEST GEOMETRY AND CONDITION INPUT PARAMETERS 


Input 

Variable 

Default 1 "' 


Name 

Type 

Value 

Description 

HTR 

R 

0.44 

Hub-tip ratio. 

OR 

R 

5.0 

Outer radius of duct. 

EMX 

R 

0.07 

Axial Mach number (always positive). 

IRQ 

R 

3)00. 

Test frequency (Heru) 

SPEED 

R 

13566. 

Speed of sound 

X(l) 

R 

9.568 

Axial coordinates of the measurement 

X(2) 


6.582 

microphone locations. (Input LOC value) 

X(3) 


0.0 

- 


X(50) 


0.0 


R(l) 

R 

5.0 

Radial coordinates of the measurement 

R(2) 


5.0 

microphone locations. (Input NLOC value) 

R(3) 


0.0 



R(50) 


0.0 


THM(l) 

R 

0.0 

Circumferential coordinates of the measurement 

THM(2) 


0.0 

microphone locations (degrees). 

THM(3) 


0.0 

(Input NLOC value) 


THM(50) 


00 

Bt 'AM(t) 

R 

0.03 1 36 Pressure amplitude at the measurement 

BE f AM(2) 


0.02097 microphone locations. (Input NLOC value) 

BETAM(3) 


0.0 

BETAM(50) 


0.0 
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TABLE II (Cont’d.) 


Input 

Variable 

Default 


Name 

Type 

Value 

Description 

PS1M(I) 

R 

97.8 

Pressure phase at the measurement micro- 

PS1M(2) 


215.6 

phone locations (degrees). 

PSIM(3) 


0.0 

(Input NLOC values) 


PS1M(50) 


Ox 


PX(!) 

R 

5.788 

Axial coordinates of the prediction micro- 

PX(2) 


2.513 

phone locations. (Input LOCP values) 

PX(3) 


0.0 



PX(50) 


0.0 


PR(I) 

R 

5.0 

Radial coordinates of the prediction micro- 

PR(2) 


5.0 

phone locations. (Input LOCP values) 

PR(3) 


0.0 



PR(50) 


0.0 


THP(I) 

R 

0.0 

Circumferential coordinates of the prediction 

THP(2) 


0.0 

microphone locations (degrees). (Input LOCP 

THP(3) 



values) 


THP(50) 


0.0 


ICHK 

I 

0 

Mode amplitude and phase indicator. 




0 - Calculated from measured pressure 




1 = Input values 

AM(1) 

R 

0.0 

Mode amplitude. (If ICHK « 1, input NLOC 

AM(2) 


0.0 

values) 

AM(3) 


0.0 


AM(50) 


0.0 
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TABLE II (Cont’d.) 


Input 

Variable 

Default 

Name 

Type 

Value 

PH1(1) 

R 

0.0 

PH1(2) 


0.0 

PHI(3) 


0.0 

PHI(50) 


0.0 


Description 


Mode phase (degrees). (If 1CHK * 1 , input 
NLOC values) 



y: 



Note, (a) Default values shown in table are in units of »he English System. The program, however, is 
designed to be executed with data in any consistent system of units. 








TABLE III 

ERROR DEVIATION COEFFICIENT INPUT PARAMETERS 


Input 

Name 

Variable 

Type 

Default 

Value 

Description 

S1GX 

R 

0.0 

Standard deviation of the axial coordinate 
error. 

S1GR 

R 

0.0 

Standard deviation of the radial coordinate 
error. 

S1GT 

R 

0.0 

Standard deviation of the circumferential 
coordinate error (degrees). 

SIGB 

R 

0.0 

Standard deviation of the pressure amplitude 
error. 

SIGP 

R 

0.0 

Standard deviation of the pressure phase 
error (degrees). 
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of degrees. The corresponding mode indices, axial wave number in units of degrees-per- 
length, and eigen value k m#i ° are delineated. 

Additional input parameters listed in this section include the reference location usually cor- 
responding to the fan face where the modal phases are calculated. Coordinates of the input 
measurement locations and resultant prediction locations are listed adjacent to the respec- 
tive acoustic pressure values. The input pressure values arc supplied by the user in pressure 
units for the amplitude and degrees for the phase. The resultant pressure is calculated by 
the program and output is provided in units of decibels for the amplitude and degrees for 
the phase. 

The Sensitivity Coefficient Calculation portion of the output comprises a number of sections, 
the primary output of which is the total normalized amplitude and the total phase deviation 
for each mode. These expressions represent the modal amplitude and phase error caused by 
a specified set of independant errors associated with the measurement of acoustic pressure 
and the tolerance of pressure measuring coordinates. The amplitude standard deviation of 
a specific mode is expressed as both the normalized quantity with respect to the mode am- 
plitude and the mode amplitude error in decibels. The total phase deviation is expressed in 
degrees for each mode. 

The contribution to the total amplitude and phase deviation assuming zero errors for the 
other error sources is provided under the heading "Normalized Standard Deviation Com- 
ponents”. The amplitude deviation was normalized with respect to the mode amplitude. 

The total phase deviation in degrees for each error source is also provided under the heading. 
When these values are root-sum-squared, the previous expression for the total modal devia- 
tion is obtained. A user will benefit from the error deviation components by identifying 
which of the errors is controlling the total modal error. 

The standard deviation components are also normalized with respect to their respective 
error. These parameters - referred to as the root-sum-square of the influence coefficients 
— enhance the combined variance of the influence coefficients at each microphone location. 
Thus, these parameters are the previous standard deviation components with respect to a 
unit measurement error in pressures or microphone coordinates. The root-sum-square of the 
influence coefficients is a convenient expression for assessing the probability of successfully 
tracking modes. A future user could examine these parameters to determine if the accuracy 
of experimental measurements made during an earlier test is sufficient to provide a desired 
confidence level in the mode amplitude and phase. 

The influence coefficients are the partial derivatives of the mode amplitude and phase with 
respect to an error at each pressure measurement location that provides input for calculating 
the modal coefficients. These expressions allow a future user to evaluate the effect of non- 
uniform errors at the microphone locations. For example, an amplitude measurement error 
may be known to be significantly larger at one microphone location (e.g.. inaccurate calibra- 
tion) The user could then evaluate the impact of this error on the overall modal structure 
calculation. 
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The final section, Characteristic E-Functions, includes the value of E-functions, E(k° r), 
at the measurement and prediction locations corresponding to each mode. This final section 
is provided as output only if it has been requested by the user. 

3.5.2 Sample Cases 

Two cases are presented in the sample printout, to illustrate the two options: I ) calculating 
mode amplitude and phase values from acoustic pressure signals and 2) specifying these 
values either arbitrarily or as output from an analytical prediction deck. These sample cases 
demonstrate the execution of each option with data listed in Figure 6. The length units in 
the printout are in centimeters; the time units, in seconds; :he force units, in dynes. 

The first sample case illustrates the option of calculating the mode amplitude and phase for 
a situation where three modes are propagating in a half-meter diameter annular duct. Three 
coherent acoustic pressure amplitude and phase values are specified at three microphone loca- 
tions on the duct wall. These acoustic signals are at a frequency of 6200-Hertz, and are used 
to determine the modal structure of the (-4,0), (-4,1) and (-4,2) modes. 






The output for this sampi case reveals that the amplitudes of the above modes are 137.4, 
142.8, and 138.5 decibels, , 'spectively; the modal phases are, respectively, 126.9, 160.0, and 
229.2 degrees. Once the modal structure has been determined, the resultant sound field can 
be calculated at other duct locations. The resultant amplitude and phase - expressed in the 
same units as the modal coefficients - are requested at three microphone coordinates. The 
resultant amplitude at these locations are, respectively, 121.1, 1 15.7, and 1 15.1 decibels. 

The resultant phases are, respectively, 90.1, 357.8, and 67.2 degrees. 

The sensitivity coefficient calculation portion of the program calculates the accuracy of the 
mode amplitude and phase values based on inaccuracies in the measured acoustic pressures 
and the microphone coordinates. Errors in the five measured quantities are expressed as 
standard deviations with zero mean. For this sample case they are axial - 2x 10*^ cm, radial - 
2x 10"’ cm, circumferential - 2x 10"^ degree, amplitude - 25 dynes, and phase - 1.5 degrees. 
The combined effects of the error source deviations multiplied by the influence coefficients 
yields the modal amplitude and phase deviation. These calculated values for the (-4, 0), (-4, 1), 
and (-4, 2) modes are, respectively, 0.89, 0.84. and 0.87 decibel for the modal amplitude 
and 3.6, 2.5, and 1 .7 degrees for the modal phase. 

The second sample case illustrates the option to input the amplitudes and phases for the pro- 
pagating modes to calculate the resultant acoustic pressure at specified locations. This case 
is similar to the first sample case because the (-4, 0), (-4, 1 ), and -4, 2) modes are propagating 
at 6200-Hertz in a half-meter diameter annular duct. The amplitude of all the modes is 
1 2 1 .9 decibels and the phases of these modes are, respectively , 325, 250, and 100 degrees. 
Output from the Modal Calculation Program comprises the resultant sound field at three 
microphone locations. The value of the resultant sound field is 1 15.8, 120.0, and 120.0 
decibels for the resultant amplitude and 36.4, 345.0, and 298.2 degrees for the resultant 
phase. 
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UTILITY CODING FORM 











3.6 MACHINE REQUIREMENTS 


T\e Modal Calculation Program can be compiled, linkage edited, and executed in 5 1 2 bytes 
of core storage. 


The following mathematical functions and procedure are required: 



CMPLX - 

Expresses two real arguments in complex form. 


CABS - 

Modulus of a complex argument. 


CEXP - 

Exponentiation of a complex argument. 

•aJW 

AIMAG - 

Obtain imaginary part of a complex argument. 


REAL - 

Obtain real part of a complex argument. 


FL0AT - 

Conversion from integer to real. 


IFIX 

Conversion from real to integer. 


ABS 

Absolute value of a real number. 


IABS 

Absolute value of an integer. 


SORT - 

Square root of a real value. 

\ 

MAXO - 

Obtain maximum value of input integers. 

1 

AL0G 

Natural iogarithim of a real positive argument. 

SIN 

Sine of a real argument. 


COS 

Cosine of a real argument. 


AT AN 2 - 

Arc tangent of two real arguments. 


3.7 RESOURCE ESTIMATES 


The central-processor-unit (CPU) time required to process a particular case depends on the 
number of modes input wnich determines the size of the matrix to be inverted. The average 
esti ..ate of CPU time per mode is 0.1S second. 


(M 
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APPENDIX A 


Calculation of the Characteristic Numbers 

The characteristic numbers IC 0 and Q£L are defined to be the paired roots of the simultaneous 
.. mu "V 1 

equations 

[A. J m (KmV) + Q r ^rY m (K° r')l r '=l *0 (1) 

dr m mu dr m mu r 1 


d * o o d 'o 

I7-7J m< oK f ) + Q T7 Y m (aK r >1 r' = i s ° 

dr m mu mu dr m mu 


For a given circumferential mode number, m, radial order, u, and hub/tip ratio, o, (where 
a is not equal to zero); J m and Y m are the Bessel functions of the first and second kinds of 
order m. 

The following relations are used in the formulation of a solution 


-filnW&Wjp- 


— Y (x)=Y M— 
df , i m w I m^ x ^ dr ' 


Vl W JmW- J m-l^) 


J m 'w-T [ V|W- J m + l Wl 


— <*> J m M * J m-1 <*>• 


■V|W T ln,W 
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. 2 
Y m (x) = 


m »xJ m (x) 


+ 4 Cx) 


Y m (x) 




2 . / , m , , Y m W 

' ,xJ m (x) 1 "- 1 <x) x m <«>' Cm 


a o . 

Letting K * K and Q = Q , and evaluating at r = 1 ; ( 1 ) and (2) become 

m/i 


y (K)K + QYL(K)K=0 
m m 


y m (oK)oK + QY^(oK)oK*0 


m 


J' m (K)K 

From (8), Q = - — t — — — — substituting into (9) yields 
Y m (K) K 


V flO K 

*m ‘oK) oK - Y' m (oK> (K=0 

Ut f (K) * l' m (oK> Y^K) v K 2 - 4 (K) (oK) oK 2 - 0 
Using the expressions in (5), (6), (7), and (11) then: 


(7) 


( 8 ) 

(9) 


(10) 

( 11 ) 


► 

f 


- t 2 U.n-1 (K) - =-i ra (K)l j -i—j- ♦ V™ J££|“0<P-> 


f(K) 


, , k 2 j 2 |J n.-l (g|Q- & V°*) l - 2|1 m-l (K> -g J m 

t — If t /V\ 1 / t / . 


(K)J 


»KJ m (K) 


xoKJ m (j K) 




L^JL’lU 

V‘ K > ( 


(13) 
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Equation (13) is evaluated for values of Kj = M + 3(i-l) ; i * 1, 2, 3, . . . until f (Kj) f (Kj * 1) < 0 
for some j. A procedure employing a combination of false position and bisection techniques 

f 

is then used to obtain a value of K 0 in the interval Kj]. 

» a 

Having calculated a value of K * K^, the corresponding value of Q = Q can be calculated. 
Combining (8) and (9) yields. 


[r m (K) + r n (oK) o J K + Q ( Y^(K) + Y' m (oK) ol K = 0 


from which 


Q» - 


J m ( K) J m (oK)a 
Y^CKl + Y^ (oK) a 


For o = 0, Q mu = 0 and K ° - 0 is defined to be the root of 


h~? J m (K° r' )1 r » = j = 0 
dr 1,1 mp 11 


Letting K = K 0 , and evaluating at x = 1,(16) becomes 
nyt 


If «K) = J ' (K) K = 0, then (6) yields 


f(K)=U m .i(K)- T J m (K)l K = 0 


Equation (18) is evaluated for values of Kj = m + 3(i-l);i = 1, 2, 3, . . . until f (Kj) f (Kj.j) < 0 
for some value of j. A procedure employing a combination ot false position and bisection 

techniques is then used to obtain a value of K ° in the interval (K : ■ , K s ). 

mu J -1 J 
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APPENDIX B 

MODAL CALCULATION PROGRAM 
PROGRAM LISTING 


blanks 


FILJdEE* 


31 


PRATT C WHITNEY AIRCRAFT DIVISION 
SC.PANLIB.L4 


VER 

9.0 


07/25/77 

12.S0.00 



♦♦WRITE PRINT* T89902 

C DATA SET T89902 AT LEVEL 024 AS OP 06/22/77 

C DATA SET T89902 AT LEVEL 014 AS OF 02/01/77 00001 

C 00002 

C THIS PROGRAM CALCULATES MOOAL AMPLITUDE ANO PHASE 00003 

C 00004 

C OOOOS 

COMMON /PREDCT/ XPI50I* RPI50I* THETAPtSOI 00006 

CCMMCN /APHIMU/ AMUI50)* PH1MUIS0I* 1CHECK 00007 

COMMON /REFCON/ REFPRS 00008 

COMMON /EMUS/ EMUl 50*501, EMUP(S0»50)* EMUPRMI50.50 I* 1ENPRT 00009 

COMMON /CNSTNT/ NMEAS* NPRED* NMUDtS, SIGMA* 8* MX* FREQ* A* 00010 

1 OMEGA 00011 

COMMON /KUMU/ KHUI50I* 0MUI50) 00012 

COMMON /MOOES/ M0DEI50), MUI50I* IWAVEISO) 00013 

COMMON /ANGLES/ OEGRAO* RAOOEG 00014 

COMMON /OUTPUT/ AMPRI50), PHASERISOl 00015 

COMMON /WAVENO/ KX ISO) 00016 

COMMON /REPS/ XREFt RREF* THREF 00017 

REAL KMU* MX 00018 

COMPLEX KX* EXPNT. SUMl* 0(50*50) 00019 

DIMENSION EMUDUMI50) 00020 

C 00021 

C INPUT DATA FOR THIS CASE 00022 

C . 00023 

20 CALL INPUT! 1EN0 ) 00024 

IFI 1END .GT. 01 GO TO 9999 00025 

C 00026 

C CALCULATE THE CHARACTERISTIC NUMBERS KMU ANO QMU FOR EACH SET OF 00027 

C CIRCUMFERENTIAL MOOE NUMBER ANO RADIAL ORDER 00029 

C 00029 

CALL KQCAL 00030 

C 00031 

C CALCULATE AXIAL WAVE NUMBER 00032 

C 00033 

FLOW « OMEGA / A 00034 

AMACH « l. - MX • MX 00035 

DO 40 I>1*NM00ES 00036 

RADI CL * FLOW ** 2 - AMACH • I KHUII) /II H 2 00037 

IFI RADICL ) 25* 30* 30 00038 

25 KXII) * CMPLXI -MX • FLOW / AMACH* lWAVE(l) • 00039 

1 SQRTI ABSI RADICL ) ) / AMACH ) 00040 

GO TO 40 00041 

30 KXII) « CMPLXI I -NX • FLOW ♦ IWAVEII) • SQRTI RAOICL ) ) /00042 

1 AMACH* 0.0 ) 00043 

40 CONTINUE 00044 

C 00045 

C IF THIS IS a CHECK RUN* AMU AND PHIMU HAVE BEEN INPUT. Thus THERE IS 00046 

C NO NEED TO CALCULATE THEM. 00047 

C 00048 

1 F ( ICHECK .GT. 0 ) GO TO 60 00049 

C 00050 


C SET UP ANO SOLVE THE EQUATION SYSTEM ASSOCIATED WITH THE MEASUREMENT 00051 
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MATT t WHITNEY AIRCRAFT DIVISION VER 07/25/77 

SC.PANLIB.L6 0.0 12.50.00 

LOCATIONS 00052 

00053 

CALL SOLVE 0003A 

00055 

CALCULATE SENSITIVITY COEFFICIENTS 00036 

0005T 

CALL SENSTVI Q* NHOOES I 00056 

00059 

CALCULATE SUN OF NODAL ANFLITUOES AND PHASES FOR EACH PREDICTION 00060 

LOCATION 00061 

00062 

60 00 120 J*1»NPRED 00063 

RPR1NE * RP( J) / 0 00066 

C 00063 

C CALCULATE CHARACTERISTIC E -FUNCTIONS FOR RPRIHE 00066 

C 0006T 

CALL ENUCALI RPRIHE* EHUPll.J). ENUOUN. 0 I 00063 

C 00069 

OXP • XPU) - XREF 00070 

OTHETP • THETAPO) - THREF 00071 

SUM1 * CNPLXI 0. ( 0. ) 00072 

DO 100 IsI.NNODES 00073 

EXPNT * CNPLXI 0.0* REAL! KXtl) I • OXP ♦ NOOEIII • OTHETP 00076 

1 • PHI NUIl) ) 00075 

SUH1 - A HU II) • EHUPll.J) • CEXPI EXPNT ) • EXPI -OXP • 00076 

1 AIMAGI KXtl ) ) ) ♦ SUM 00077 

100 CONTINUE 00076 

C 00079 

AHPRIJ) - CA6SI SUNl ) 00080 

PHASERIJI * ATAN2I AINA6I SUNl ) ( REALI SUNl t I 00061 

120 CONTINUE 00082 

C 00083 

C PRINT RESULTS OF THIS CASE . 00066 

C 00065 

CALL PRINT 00096 

C 00067 

C RECYCLE FOR NEXT CASE 00088 

C 00069 

60 TO 20 00090 

9999 STOP 00091 

END 00092 

SUBROUTINE ANGPOSI ANGLE. NUHOER ) 00093 

C 00096 

C THIS SUBROUTINE CONVERTS NEGATIVE ANGLES TO CORRESPONDING POSITIVE 00095 

C ANGLES 00096 

C 00097 

DINENSION ANGLE! I) 00098 

DATA DEGREE / 360. / 00099 

C 00100 

DO 60 1*1*NUNBER 00101 

1F( ANGLE! II ) 20* 60* 60 00102 

20 DO 60 J*1.10 00103 

DELTA ■ J • DEGREE 00106 
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PRATT & WHITNEY AIRCRAFT DIVISION 
SC.PANLIb.L4 


VER 

9.0 


07/25/77 

12.50.00 




IF( ANCLE I II ♦ DELTA I 40* 60. *0 00105 

40 CONTINUE 00106 

60 ANGLel 1} > DELTA ♦ ANCLE! I I 00107 

60 CONTINUE 00108 

9999 RETURN 00109 

6 Nil 00110 

FUNCTION BESLli X I 00111 

C 00112 

C THIS FUNCTION CALCULATES VALUES OF THE EQUATION DEFINING THE SYSTEN 0F00U3 

C DIFFERENTIAL EQUATIONS FOR A NON-ZERO HUB/TIP RATIO 00114 

C 00115 

COMM IN /BESSL/ ISIGN» JSIGN. OELKNUt TOL. N» PI 00116 

COMMON /CNSTNT/ DUM1I3). SIGMA. 0UM215I 00117 

C 00116 

XI ■ X • SIGMA 00119 

CALL BESJ I XI. M-JSIGN. EMJM1, TOL, IER1 I 00120 

CALL BESJI XI* N * EMJXl. TOL. 1ER2 I 00121 

CALL BESJI X. M. EMJ. TOL. IER3 I 00122 

call besji x, m-jsicn. emjpi. tol* iera i 00123 

CALL BESYi X. M, EMYX. 1ER5 1 00124 

CALL BESYI XI* M. EMVX1. 1ER6 I 00125 

C 00126 

EMJH1 * JSIGN • 1SIGN * EMJN1 00127 

EMJX1 ■ I SIGN • EMJXl 00126 

EMJ * I SIGN • EMJ 00129 

EMJPI b JSIGN • I SIGN • EMJPI 00130 

EMYX - 1SIGN * EMYX 00131 

EMYX 1 - 1SIGN • EMYX1 00132 

C 00133 

A1 ■ EMJM1 -IN* JSIGN / XI » * EMJXl 00134 

A2 « EMJPI -IN* JSIGN /X I • EMJ 00135 

A3 * 2. * A1 / I PI • X • EMJ I 00136 

A4 « 2. • A2 / I PI • XI • EMJXl 1 * 00137 

A5 ■ Al • A2 • I EMYX / EMJ - ENVX1 / ENJX1 I 0013B 

C 00139 

BESL1 » XI * X • I A3 - A4 • A5 ) 00140 

RETURN 00141 

END 00142 

FUNCTION BESL2I X ) 00143 

C 00144 

C THIS FUNCTION CALCULATES VALUES OF THE EQUATION DEFINING THE SYSTEM 0F00145 

C DIFFERENTIAL EQUATIONS FOR A HUB/TIP RATIO OF ZERO 00146 

C * 00147 

COMMON /BESSL/ IS1GN. JSIGN. OELKHU. TOL. M. PI 00148 

COMMON /CNSTNT/ DUMIOI. SIGMA. 0UM2I5I 00149 

C 0015P 

CALL BESJI X. M-JSIGN. EMJMl . TOL* lERl I 00151 

CALL BESJI X. M, EMJ. TOL. IER2 ) 00152 

C 00153 

EMJMl > JSIGN • I SIGN • EMJMl 00154 

EMJ * 1SIGN * EMJ 00155 

BESL2 b x • EMJMl - N • JSIGN • ENJ 00156 

RETURN 00157 
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PRATT C WHITNfcV AIRCRAFT DIVISION VCR 07/25/77 *AGE SERIAL 

SC.PANL1B.L4 9.0 12.50.60 ** 010B76 


f NO 00159 

SUBROUTINE BESJI X, N. 8J. Ot IER I 00159 

C OOlbO 

C THIS SUBROUTINE CALCULATES THE J BESSEL FUNCTION FOR A GIVEN ARGUHENT»00161 

C X* ANO ORDER N. THIS SUBROUTINE NAS TAKEN FROM THE IBM SCIENTIFIC 00162 
C SUBROUTINE PACKAGE 00163 

C 00164 

BJ ■ 0.0 00165 

1FI N .GE. 01 GO TO 20 00166 

C 00167 

C ERROR - NEGATIVE ORDER. SET ERROR 1N01CAT0R TO 1 ANO RETURN 0016B 

C 00169 

IER * 1 00170 

GO TO 9999 00171 

20 IFI X ) 40. 30. 60 00172 

30 IFI N .GT. 01 GO TO 40 00173 

BJ * 1.0 00174 

GO TO 9999 00175 

C 00176 

C ERROR - ARGUMENT ZERO OF NEGATIVE. SET ERROR INOICATOR TO 2 AND RETURN00177 
C 00178 

40 IER * 2 00179 

GO TO 9999 001B0 

C 001B1 

C CALCULATE MAXIMUM ORDER NUMBER THAT CAN BE PROCESSED FOR X. 00162 

C IF X .LE. 15. N MUST BE LESS THAN 20 ♦ 10*X - X**2/3. 001B3 

C IF X .GT. 15. N MUST BE LESS THAN 90 ♦ X/2 00164 

C 00185 

60 IFI X - 15. I 60. 80. 100 00186 

90 NTEST * 20. ♦ 10. • X - X •* 2 / 3. 00187 

GO TO 120 00199 

100 NTEST « 90. ♦ X / 2. 001B9 

120 IFI N .LT. NTEST 1 60 TO 140 * 00190 

C 00191 

C ERROR > ORDER RANGE COMPARED TO X IS NOT CORRECT. SET ERROR INOICATOR 00192 

C TO 4 AND RETURN. 00193 

C 00194 

IER * 4 00195 

GO TO 9999 00196 

140 IER ■ 0 00197 

N1 * N ♦ 1 00198 

BPREV » 0.0 00199 

C 00200 

C COMPUTE STARTING VALUE OF N 00203 

C 00202 

IFI X - 5. I 160. ISO. 180 00203 

. 160 MA ■ X ♦ 6. 00204 

GO TO 200 00205 

180 MA * 1.4 • X ♦ 60. / X 00206 

200 MB « N ♦ 1FIXI X » / 4 ♦ 2 00207 

MZERO > MAXOI MA. MB I 00208 

C 00209 

C SET UPPER LIMIT OF M 00210 
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PRATT & WHITNEY AIRCRAFT OlVISION 
SC.PANUft.L4 


VER 

9.0 


07/25/77 

12.50.00 



C SET F !M ) • FIM-D 
C 

FMl 

FM 

ALPHA 

JT 


1.0E-25 

0.0 

0.0 

1 


IF! ! H / 2 I • 2 *EQ. H I 
M2 « M - 2 

DO 280 Mi#M2 
HK « M - K 


JT * -I 


IF! MR - N - I I 
bj * bMK 


240 bJ 
260 JT 
S 

ALPHA 

260 CONTINUE 


M R 

2. ft FLOAT! MR I ft FMl / X - FN 

f MI 

ftfftC 

) 260* 240* 260 


C 00211 

MMAX ■ NTEST 00212 

220 00 320 MbMZER0.MMAX»3 00213 

C 00216 

C SET FINI, FIN«1) 0021S 

C 00216 

FMl - 1.0E-2B 00217 

FM > 0.0 0021U 

ALPHA • 0.0 00219 

JT >1 00220 

IFI < M / 2 > • 2 .EQ. M ) JT * -1 00221 

M2 ■ M - 2 00222 

DO 260 K*1»M2 00223 

MK » M - K 00226 

BMK - 2. • FLOAT! MK ) • FHl / X - FN 00226 

FM « FMl 00226 

FMl > BMK 00227 

1FI MK - N - 1 ) 260* 260* 260 00228 

260 BJ * bMK 00229 

260 JT « -JT 00230 

S ■ 1 * JT 00231 

ALPHA - ALPHA ♦ BMK • S 00232 

280 CONTINUE 00233 

C , 00236 

BMK ■ 2. • FMl / X - FM 00236 

IFI N .tO. 01 BJ • BMK 00236 

ALPHA • ALPHA ♦ BMK 00237 

BJ > IJ / ALPHA 00238 

IF! ABSI BJ - BPREV I - ABSI 0 * BJ I I 9999* 9999* 300 00239 

300 BPKtV ■ BJ 00260 

320 CONTINUE 00261 

C . 00262 

C ERROR - REQUIRED TOLERANCE NOT OBTAINED. SET ERROR 1N01CAT0R TO *3 AND 00263 
C RETURN 00266 

C 00265 

1ER - 3 00296 

9999 RETURN 00267 

END 00269 

SUBROUTINE BEST! X* N* BY* IER ) 00269 

C 00260 

C THIS SUBROUTINE CALCULATES THE V BESSEL FUNCTION FOR A GIVEN ARGUMENT *002 51 

C X* ANO UtDtR N. THIS SUBROUTINE MAS TAKEN FROM THE IBM SCIENTIFIC 00252 


-JT 

1 * JT 

ALPHA ♦ BMK • S 


BMK ■ 2. • FMl / X - FM 

IFI N .EO. 01 BJ • BMK 

ALPHA • ALPHA ♦ BMK 

BJ * BJ ✓ ALPHA 

IF! ABSI BJ - BPREV I - ABSI D * BJ I I 9999* 9999* 300 
300 BPREV b BJ 

320 CONTINUE 


SUBROUTINE PACKAGE 

IER b C 

IFI N .GE. 0 ) 


GO Tq 20 


C ERROR - NEGATIVE ORDER. SET ERROR IN01CAT0R TO 1 ANO RETURN 
C 

IER b 1 


20 IFI X ) 


GO TO 9999 
60* 60* 60 


00253 

00266 

00265 

0026b 

00257 

0026B 

00269 

00260 

00261 

00262 

00263 
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SC.PAMLIS.L4 


VER 

*•0 


07/25/77 

I2.S0.00 


EAROft - ARGUMENT ZERO OR NEGATIVE. SET ERROR INDICATOR TO 2 AND RETURN00244 

0024 S 

40 1ER *2 

GO TO VOW 


C SRANCH IF X IS LESS THAN OR EQUAL TO 4. 
C 


40 1FI X - 4. I 


100. 100, 00 


C CALCULATE VO AND VI FOR X GREATER THAN 4. 

C 

- 4. / X 
• T1 # T1 

■ I I I I -.0000037043 #TI ♦ .00001TSS4S » # T2 - 00277 

• 0000407413 I #T2 ♦ .00017343 1 • T2 - .001753042 0027S 


SO T1 
T2 
PO 

1 

2 

00 

1 

2 

PI 

1 

2 

01 

1 

2 

A 

B 

C 

VO 

VI 


00244 
00247 

00245 
0024* 

00270 

00271 

00272 

00273 

00274 

00275 
00274 


• • T2 ♦ .3*0*423 

■ I ( < I .0000032312 # T2 - .0000142070 > * T2 ♦ 
•0000342440 » # T2 - .0000*4*7*1 I # T2 ♦ 
.0004544324 I # T2 - .012444*4 

■ I I I I .0000042414 # T2 - .0000200*2 I # T2 * 
.0000500759 I • 72 - .000223203 | # T2 ♦ 
.002921024 I • T2 ♦ .3*0*423 

> 1 1 1 1 -.0000034594 * T2 ♦ .00001422 I ATI • 
.0000390700 I # T2 .0001044741 » # 72 - 
•00043904 I #T2 ♦ .03740004 

> 2. / SORT I X I 

• A • T1 

> X - .7053902 

• A • PO ♦ SINI C ) ♦ 0 • 00 * COSI C ) 

> -A • PI # COSI C I ♦ 0 #01 AtlNI c > 

GO TO 140 


C CALCULATE VO AND VI FOR X LESS THAN OR EOUAL TO 4 * 

C 

100 XX > .5 • X 

X2 > XX • XX 

T »< ALOGI XX I ♦ .5772157 

SUN * 0.0 

TERN * T 

VO ■ T 

00 120 L>1»15 

IF< L .HE. 11 SUN > 1. / FLOAT I L-ll ♦ SUN 

FL > L 

TS > T - SUN 

TERN • I -X2 • TERN / I FL •• 2 I ) • I 1. - 1. ✓ I FL • 

1 TS I I 

VO * TERN ♦ VO 

120 CONTINUE 

TERM > XX • | T - .5 I 

SUN * 0.0 

VI * TERN 

00 140 L>2«14 

SUN • 1. / FLOAT* l-l I ♦ SCO* 


00279 

00200 

00201 

00202 

00263 

00204 

00205 
00204 
002*7 
00200 
0020 * 

00290 

00291 

00292 
002*3 
00294 

00296 
00294 

00297 
00290 
0029* 

00300 

00301 

00302 

00303 

00304 

00305 
00304 
00307 
00300 
0030* 
00310 

00311 

00312 

00313 

00314 

00315 
00314 
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PRATT & WHITNEY AIRCRAFT DIVISION 
SC*PANLIB.L4 


VI 

140 CONTINUE 
PI2 
VO 
VI 


■ L 

■ PL - 1. 

* T - SUN 

* I -*2 • TEAM / I El • FU M • I I TS - .5 ✓ PL I 
/ t T$ • .5 / PLl I I 

• TERN ♦ VI 

• .6)66190 

■ PI2 ♦ VO 1 

■ PI2 • I Y1 - 1* / I | 


C CHECK IF ONLY VO OR VI IS DESIRfO 
C 

160 1FI N «GT. 1 | 

C 

C RETURN VO OR VI AS REOUIREO 
C 

6V • VO 

FI N .EC. 1 I 


GO TO 100 


0V • VI 
GO TO 9999 


C PERFORN RECURRENCE OPERATIONS TO FINO YN(X) 

C 

100 VA « VO 

Y0 « VI 

K « 1 

200 T * FLOAT I ?*K I / X 

VC * T * V0 - VA 

IF! ABSI VC > - 1.0E7Q I 240, 240, 220 

ERROR -BY HAS EXCEEDED MACN1TU0E OF 109970. SET ERROR INDICATOR TO 3 
AND RETURN 

f 

220 1ER ■ 3 

GO TO 9999 

240 K « 1 ♦ K 

IF I K «EQ. N I GO TO 260 

VA « Y0 

YB » VC 

GO TO 200 

260 BY « VC 

9999 RETURN 
END 

BLOCK DATA 

COMMON /OFAULT/ LOCM, LOCP, HTR, OR, ENX, FRQ, XO, RO, 

1 THO, XI501, RI50), THMI50), BETAM(50), PSIHI50), 

2 PX150I, PR ( 501 1 THPOO), M<S0I, NUS(50), 

3 IDIR(SO), PREF, ANISO), PHI f 50 It ICHK, 

4 SIGX , S1GR , SIGT, SIG0, SIGP, lEMUt SPEEO 

CONSTANT DEFAULT VALUES 

LOCN - NUMBER OF MEASUREMENT LOCATIONS 


00317 

00310 

00319 

00320 

00321 

00322 

00323 

00324 
003J:6 

00326 

00327 
00329 

00329 

00330 

00331 

00332 

00333 

00334 

00335 

00336 

00337 
00336 

00339 

00340 

00341 
00362 

00343 

00344 

00345 

00346 

00347 
00340 

00349 

00350 

00351 

00352 

00353 

00354 

00355 

00356 

00357 

00358 

00359 

00360 

00361 

00362 

00363 

00364 

00365 

00366 

00367 

00368 

00369 
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C LOOP - NUMBER OP PREDICTION LOCATIONS 00370 

C NTS - HUB / TIP HAT10 00371 

C OK OUTER RADIUS 0037? 

C ENX - AXIAL MACH NUMBER 00373 

C PRO - TEST FREQUENCY 00374 

C SPEED - SPEED OF SOUNO 00373 

C 00376 

OATA LOOM / 2 /, LUCP / 2 /• HTR ✓ 0.46 /, OR / 5.0 /, 00377 

1 ENX / -*.07 /, FRO / 3100. /. SPEED / 13566.24 / 00378 

00379 

REFERENCE LOCATION VALUES ' 003UU 

00381 

XO - AXIAL CUMPuNENI OF REFERENCE LOCATION 00382 

HO - RAOUL COMPONENT OP REFERENCE LOCATION 00383 

7 NO - ANGULAR COMPONENT OF REFERENCE LOCATION (0E6I 00364 

00383 

OATA XO / 0.0 /, RO / 0.0 /, THO / 0.0 / 00386 

00367 

MEASUREMENT LOCATION VALUES 00368 

00389 

X - AXIAL COMPONENT OF MEASUREMENT LOCATION 00390 

R - RADIAL COMPONENT OF MEASUREMENT LOCATION 00391 

TON - ANGUIAR COMPONENT OF MEASUREMENT LOCATION I0E6I 0039.' 

BETAM - AMPLITUDE OF MEASURED VALUE , 00393 

PSIM - PHASE ANGLE OF MEASURED VALUE IOEOI 00394 

00393 

DATA X / 9.368, 6.382, 48*0.0 /, R / 2*3.0, 48*0.0 /, 00396 

1 THM / 30*0.0 /, bfcTAM / 0.63136, 0.05097, 48*0.0 /, 00397 

2 PSIM / 97.8, 213.6, 48*0.0 / 0C398 

00399 

PREDICTION LOCATION VALUES Of 400 

06*01 

PX - AXIAL COMPONENT OF PREDICTION LOCATION 0040? 

PR - RAOUL COMPONENT OF PREDICTION LOCATION 00403 

THP - ANGULAR COMPONENT OF PREDICTION LOCATION IOEOI 00404 

00403 

DATA PX / 3.788, 2.313, 48*0.0 /, PR / 2*5.0, 48*0.0 /» 00406 

1 THP / 30*0.0 / 00407 

00408 

MOOE VALUES 00409 

00410 

N - CIRCUMFERENTIAL MOOE NUMBER 00411 

MUS - RADIAL ORDER 00412 

101R - INCIDENT CR REFLECTED WAVE INDICATOR 00413 

00414 

OATA N / >2, -2, 48*0 /, MUS / 0, I* 48*0 /, IOIR / 1* I, 48*0 / 00413 

00416 

REFERENCE CONSTANTS 00417 

00418 

PREF - REFERENCE PRESSURE 004 19 

00420 

OATA PREF / 2.9E-9 / 00421 

C 00' 2? 
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VCR 

9.0 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

t 


RESULTANT AMPLITUDE ANO PHASE VALUES 

AH - AMPLITUDE 

PHI - PHASE ANGLE (DEG) 

ichk - Check case indicator 
data AM / 50*0.0 /• PHI / 50*0.0 /, ICHK / 0 / 
BESSEL FUNCTION VALUES 

COMMON /BE SSL/ 0UM2I2)» DELKMU, TOL# MM# FI 
UAlA DELKMU / 3. ✓ * TOL / .0001 /# FI / 3.141593 / 

ANGULAR CONVERSION FACTORS 

DEGRAD - DEGREES TO RADIANS 
RAUOEG * RADIANS TO DEGREES 

COMMON /ANGLES/ DEGRAD* R A DOE 6 

DATA GEGRAD / .0174533 /# RADDEG / 57.29578 / 

ERROR DEVIATION COEFFICIENTS 


C 

C 


SIGX - AXIAL COEFFICIENT 

SIGR - RADIAL CUEFF IClLNT 

SILT - ANGULAR COEFFICIENT 

SIGB - AMPLITUDE COEFFICIENT 

S1GP - PMASL COEFFICIENT 

DATA SIGX / 0.0 /» SIGR / 0.0 /# SIGT / 0.0 /# SIGP / 0.0 /# 

I SIGB / 0.0 / 

END * 

SUBROUTINE EMUCAL I RPAIME* EMU* EH UP RM# IDERIV I 

THIS SUBROUTINE CALCULATES NHODES CHARACTERISTIC E-FUNCTION VALUES 
a particular radial value* rprime. 

DIMENSION EMUlil* EMUPRMI1) 

COMMON /CNSTNT/ NMcAS * NPREO* NHODES* 0UH1(6» 

COMMON /KOMU/ KMUI50I* QMUC50I 
COMMON /MOOES/ MODE I 50) » UUM2I100) 

COMMON /BE SSL/ ISIGN, jSiGN* DELKMU* TOL# M# PI 
REAL KMU* JPKIME 


DO 40 1 K I* NMDDES 


IABSt MOOE(I) ) 


1M M .NE. 0 I 


ISIGN 

JSlGN 

10 ISIGN 


* l 

* -I 


MOOEIII / M 


GO TO 10 
GO TO 20 


07/25/77 

12.50.00 

00423 

00424 

00425 

00426 

00427 
0G42H 

00429 

00430 

00431 

00432 

00433 

00434 
C0435 

00436 

00437 
0043 a 

00439 

00440 

00441 

00442 

00443 

00444 

00445 

00446 
C044 7 
00440 

00449 

00450 

00451 

00452 

00453 
004 54 

00455 

00456 

00457 
FOR0045B 

00459 

00460 

00461 

00462 

00463 

00464 

00465 

00466 

00467 
00460 

00469 

00470 

00471 

00472 
004/3 

00474 

00475 


t 1 


40 


i 


i 


4 i 


l 
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9.0 
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I2.S0.00 


JSIGN • ISIGN 00476 

IF I ISIGN .GE. 0 } GO 10 20 00477 

C 00470 

C NEGATIVE NOOE NUMBER. IF EVEN* SIGN OF BESSEL FUNCTION WILL BE *U IF 00479 
C ODD, SIGN OF BESSEL FUNCTION WILL BE >1. 00460 

C 004# 1 

IF( I N / 2 I • 2 .EQ. N I ISIGN » 1 00462 

20 CONST ■ KMUII) • RPRINE 00483 

C 00484 

c calculate bessel functions of first and second kind for knuii>*rprine 0046s 

C 00486 

CALL BESJI CONST* M« EMJ * TOL, 1ERI I 00467 

CALL BESVI CONST, M, ENV, 1ER2 I OO46B 

CALL BESJI CONST. M-JSIGN* ENNI* TOL* 1ER3 ) 00489 

EMMI s ISIGN • JSIGN • EMM I 00490 

EMJ « ISIGN • EMJ 00491 

E MV * ISIGN * EMV C0492 

C 00493 

C CALCULATE CHARACTERISTIC E-f UNCTION 00.94 

C 0049 S 

EMUl 1) ■ EMJ ♦ OMUII) • EMV 00496 

IFI 10ERIV .LE. 0 ) GO TO 40 00497 

1FI KMUII) > 30, 25, 30 00498 

C , 00499 

C 10,01 CASE. SET DERIVATIVE TO 0.0 00560 

C 00S01 

25 l HUP KM U» « 0.0 00502 

GO TO 40 00503 

C 00504 

C CALCULATE DERIVATIVE OF CHARACTERISTIC E~fUNCTION 0050a 

C 00506 

30 JPRIME « EMMI - MODE II) • EMJ / CONST 00507 

V PRIME * 2. / I FI • CONST • EMJ I ♦ JPRIME • ENV / Eli) 0050B 

EMUPKMI 1 1 » KMUII) • I JPRIME ♦ UMUI1) * VPRIME I 00509 

40 CONTINUE 00510 

9999 RETUHN 00511 

END 00512 

FUNCTION FAL21P IFUNCT* AL , BR , TOL* ROOT, ITER* W) 00513 

C 00514 

C CORRESPONDS TO OLD VERSION I FALSIE ) ARGUMENT LIST AS FOLLOWS ITHIS IS 00515 

C FOR INTERNAL PURPOSES ONLY, IN USE THE TWO ARE INTERCHANGEABLE). 00516 

C 00517 

C FUNCTION FALSIE IAXR, XXL* XXR, TOL, ROOT* ITER* VV) 00518 

C 00519 

C 00520 

C THIS ROUTINE USES A COMBINATION OF FALSE POSITION ANO BISECTION 00521 

C TEUiNluUES TO SOLVE FUR A ROOT l*ROOT*) OF A GIVEN FUNCTION 00522 

C I *FUNCT • ) WHICH HAS ONE ARGUMENT ITHE INDEPENDENT VARIABLE). 00523 

C 00524 

C •AL.bk* DEFINES THE INTERVAL TO BE SEARCHED. 00525 

C 00526 

C THE VALl'*. RETURNED BV THE FUNCTION IS FALI1P. FUNCTIFALZIP ) * ROOT 00527 

C 00528 
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VER 

*.0 


07 / 25/77 

I2.SO.OO 


C THt SI ARCH CONTINUES UNTIL TWO SUBSEQUENT GUESSES ARE WITHIN *T0t* 00*29 

C OF EACH OTHER. OR UNTIL 'ITER* ITERATIONS HAVE TAKEN PLACE. 00*30 

C 00*31 

C •YY* IS RETURNED AS FUNCTIFALZIP) * AND SHOULD BE CLOSE TO "ROOT*. 00*32 
C 00*33 

C THE TECHNIQUE WAS ADAPTED FROM AN ALGOL SUBROUTINE APPEARING IN THE 00*34 
COMPUTER JOURNAL 12 11969) — 'EIGENVALUES OF A*X > LAMBOA»B*X 00*3* 

W11H t AND SVHHLTR1C A AND B* BY G. PETERS ♦ J.H. WILKINSON 00*36 

00*37 

EXTERNAL FUNCT 00*3B 

REAL 1NIEKP 00*39 

00*40 

J IS COUNT OF ITERATIONS. 00*41 

1 J * 0 00*42 

A * AC 00*43 

b * bK 00*44 

00*4* 

EVALUATE FUNCTION AT LEFT (A) AND RIGHT <Bt BRACKETS. 00*46 

AF * FUNCT (A) 00*47 

BF » FUNCT (B) 00*48 

00*49 

TFiE f ULLUM1NG I THROUGH STATEMENT 31 DETERMINES IF THE FUNCTION IS OF 00**0 
UPPUS1H SIGN AT THE ENDPOINTS GIVEN. 00**1 

ISM - 1 , 00**2 

IF ( OF - ROOT) 2. 7*. 3 00**3 

2 1SW * -1 00554 

3 IF I IAF - ROOT) * 1SW) 50. SO* S* 00*5* 

00*56 

STATEMENT * INCREMENTS THE COUNTER Jt FIRST TINE THROUGH GO TO SO. 00**7 
* J » J ♦ 1 00**8 

00**9 

IF LEFT BRACKET HAS 'SAME' FUNCTION VALUE AS RIGHT. USE BISECTION. 00*60 
OTHERWISE, SET UP INTERPOLATED POINT FOR POSSIBLE USE. ' 00*61 

IF (AbSMAF - OF )/BF ) - l.E-5) 10. 10* 15 00*62 

10 1NTERP «• BISECT 00*63 

GO TO 20 00*64 

1* 1NTERP * |A*BF - B*AF ♦ IB-A)»ROOT> / IBF-AF) 00*6* 

00*66 

IF WITHIN A TOLERANCE OF THE BRACKET B. MOVE THE INTERPOLATED POINT 00*67 
ONE TOLERANCE AMAV. 00*68 

20 IF ( IAbSllNTEKP~b)/ABS(lNTERP»B)> -2.*T0L) 22,23.23 00*69 

22 INTEkP > M 1C - Bl / ABS «C - B) • TOL 00*70 

00*71 

SET A-b IB IS ALWAYS THE POINT WITH SMALLEST I ABS) VALUE OF FUNCTI0N.00S72 

23 A . B . 00*73 

AF > BF 00*74 

00*7* 

USE POINT CLOSEST TO B IINTERP OR BISECT) AS NEW B AND EVALUATE BF. 00*76 
IF II1NTEHP - BISECT) • IB - INTERP)) 30, 25, 2* 00*77 

2* B * INTERP 00*78 

GO TU 35 00*79 

30 b • BISECT 00*80 

3* bF = FUNCT IB ) G0*81 


i 


i 
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oT/as/rr 
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BFHR ■ BF - ROOT 00582 

00583 

IF CF IS UN THE SAME SIDE OF THE ROOT AS BF, LET POINT C ■ POINT A. 00584 
40 IF I (CF - ROOT I • BFMRI 55. 75. 50 00565 

50 C « A 00586 

CF =■ AF 00587 

00568 

IF CF IS CLOSER IABSI TO ROOT THAN BF* SNITCH POINTS • A NO C. 0058V 

in any case, b ano c are the tmo brackets, also BF is closer to theoosvo 

ROOT Than CF IS. 00591 

55 IF (AttSIBF - ROOT I - ABSICF - ROOTII 60. 60. 57 005V2 

57 A -to 00593 

AF » bF 00594 

B « C 00595 

6F * CF 00546 

C ® A 005V7 

CF * AF 00595 

00549 

SET UP BISECTION POINT. IF CLOSE ENOUGH. FINISH UP. 0THERM1SE GO 00600 
BACK IF ITERATION COUNT DOESN'T EXCEED MAXIMUM. 00601 

60 BISECT » IB ♦ C» / 2. 00602 

IF I TABS (BISECT— bl/ABS IBISECT+B) I -2.*T0L I 75.65.65 00603 

65 IF IJ - ITER) 5. 70. 70 00604 

70 MRlTt(6.1000)J.b.BF.L.CF , 00605 

1000 FUHMAT OHO/// 30X, 'IN FALZIP, AFTER*. 14. • ITERATIONS* // 00606 

1 101. ’bRACKEI 1>*. 015.6, 5X, 'FUNCTION a *. G15.8/ 00607 

2 luX. ’BRACKET 2 * *, 015.6* 5X, ’FUNCTION ■ *. G1S.6/ 00608 

3 5X« ’BRACKET 1 MAS RETURNED AS RESULT.’ I 00609 

75 FALZIP « B 00610 

VV » bF 00611 

RETURN 00612 

80 FALZIP A A 00613 

YV « AF ' 00614 

RETURN C"615 

85 MRI Tt 1 6. 1 100) ROOT, A. AF .B.BF 00616 

1100 FORMAT |'0***IN FALZIP. ROOT G1VFN I**, G15.8, *1 OION**T FALL BET00617 
IRE EN VALUES UF FUNCTION AT BRACKETS GIVEN****/ 00618 

2 10X, 'BRACKET 1 > *, G15.8, 5X, 'FUNCTION « •• G15.8 / 00619 

3 10X. 'bRACKLT 2 * * . G15.B, 5X, 'FUNCTION « •* G15.8 / 00620 

4 40X . 'TERMINATING RUN* I 00621 

STOP 00622 

END 00623 

SUbREJUllNE INPUT I (END I 00624 

00625 

THIS SUBROUTINE INPUTS THE DATA REQUIRED FOR THE EXECUTION OF A CASE 00626 

00627 

COMMON /DFAUL7/ COCM, LOCP. HTR. OR. EMX. FRO, XO, RO, 00628 

1 THO. X 1 501 • RI50I, THMI50I, BETAMT 50) , PS1MT50I. 0062V 

2 PXI50T, PR (50) . THPI50T, MI50), MUSI50), 00630 

3 1DIR ( 50) , PREF, AMI50), PHII50), 1CHK, 00631 

4 SlGX , S1GR , S1GT, SIGU. SIGP, 1EMU, SPEED 00632 

COMMON /CNSTNT/ NMt-AS, NPRtO, NMOOtS. SIGMA, B. MX. FREQ. A, 00633 

1 EjMEuA 00634 
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PRATT 4 WHITNEY AIRCRAFT DIVISION VCR 

SC*PANL1B*L4 M 

C OHM ON /REES/ XREF* RREF* THREF 

COMMON /NEASUR/ XHISOI* RMI50I* THETAHtSOI* liTMSOli P$ll*0| 

COMMON /PREllCT/ XPI50I* RPI50I* THETAFf SOI 

COMMON /MUOtS/ MOOt I SOI * MU I SOI. INAVEISOI 

COMMON /bt SSL/ CMJMl IS I. PI 

COMMON /RttCON/ REPPRS 

COMMON /ANCLES/ lit GRAD* RADDEG 

COMMON /APH1MU/ AHUI50I* PH1MUISOI* ICHECK 

COMMON /EMoS/ EMUl SO .SOI * EMUP (50*501* EMUPRNI50*50l* IENPRT 
COMMON /0V1AIE/ SIGH AX* Si OMAR* SIGMAT* S1GMAB* S1GMAP* DUM2I150I 
I IOEV 

C 

NAME El St /INDAtA/ LULM. EOCP. HTR* OR* EHX* FRQ* XO* 

1 TWO* X* R* THM. BEIAM* PSIM. PX. PR* THP* N* 

2 MUS* 1DIR* PREF* VREF* AM* PHI* ICHR* 

3 S1GX* SlGR* SIGT* SIGB* SIGP* I EMU. SPEED 
RE At MX* LREF 


I END » 0 

READ 1 5 .INOATA#EN0*9996 1 
C 

C SET UP INTERNAL PARAMETERS 
C 


NMEAS « 

APR tD * 

NMOOtS * 

SIGMA « 

MX * 

PktQ « 

REPPkS « 

1 CHECK « 

XREE « 

kREP • 

THREF * 

B * 

A * 

IENPRT 

SICMAX « 

SIGMAR * 

SIGMAT * 

SICMAB 
SIGMAP 

DO 20 I«i*NMEA 

xmiii 

RMill » 

THE! AMI I | « 

PS! Ill * 

HtfAllI « 

20 CONTINUE 


LOOM 

EOCP 

LOCM 

MIR 

EHX 

PRO 

PREF 

1CHK 

XO 

RG 

THO * OEGRAO 

OR 

SPEED 

I EMU 

S1GX 

SlGR 

SIGT 

SIGB 

SIGP 


X f 1 1 
Rill 

THM! I! V OEGRAD 
PSIMlIl • OEGRAO 
BETAMIII 


DO 40 I*1»NPREU 


07/21/71 

12*50*00 

0063 5 
00636 
0063T 
006 3B 

00639 

00640 

00641 

00642 

00643 
1*00644 

00643 

00646 

00647 
00649 

00649 

00650 

00651 

00652 

00653 

00654 

00655 

00656 
0065? 
0065B 

00659 

00660 
00661 
00662 

00663 

00664 

00665 

00666 
G0667 
00668 

00669 

00670 

00671 

00672 

00673 

00674 

00675 

00676 

00677 
00679 
00679 
006B0 
00661 
00662 
0066 3 
00694 

00665 

00666 
00667 
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PRATT t WHITNEY AIRCRAFT 01 VI SION 
SC.PANtl6.L4 


VIA 

9.0 


0T/25/T7 
12. SO .00 


XPI1I 
RPI1I 
THETAPII) 
40 CONTINUE 


pxm 

PRIll 

thpid 


DEGRAD 


00 40 1-ltNNOOIS 


MODE (1 1 
HUIII 
1WAVEI1) 
AMU 1 1 I 
PHlMUtll 
40 CONTINUE 


MU) 

HUSH) 

101RI1) 

AMI 1 ) 

PH! II I • DEGRAD 


DO 80 J*i,20 
DO 80 1*1,20 
ENUIltJ) * 0.0 

EMUP(l.J) * 0.0 
ENUPRHl 1 ,J) • 0.0 

80 CONTINUE 

CALCULATE KADI AN FREQUENCY 


OMEGA 


• 2. • PI * FREQ 


set indicator for error source standard deviations 


1FI SIGN AX I 
100 1FI SIGMAR ) 
120 IF I SIGMAT I 
140 1FI S1GMAB ) 
160 1FI SIGMAP I 
140 lDtV 

200 I DEV 


END OF DATA SET 


200 , 100 , 200 
200 , 120 , 200 
200, 140, 200 
200, 140, 200 
200, 180, 200 

GO TO 9999 

GO TO 9999 


9998 1EN0 • 1 

9999 RETURN 
END 

SUBROUTINE KMUCALI VALUE, DELTA, KMU, RIGHT ) 

THIS SUBROUTINE CALCULATES THE CHARACTERISTIC NUMBER, KMU 
EXTERNAL BESL1, BESL2 

common /cnstnt/ dumiui, sigma, oun2I5) 

REAL KMU, LEFT 
SO I PLUS * 0 

I MINUS « 0 

34 IFI SIGMA I 50, 40, 50 


ORIGINAL RAGBJB 
OP POOR QUALTTb 


0060 b 

00689 

00690 

00691 

00692 

00693 
00696 
00696 
00696 
0069? 

00698 

00699 
00?00 

00701 

00702 

00703 
0G7o4 
00706 

00706 

00707 

00708 

00709 

00710 

00711 
0071* 

00713 

00714 
00716 

00716 

00717 

00718 

00719 

00720 

00721 

00722 
0072 3 
00724 
00726 

00726 

00727 

00728 

00729 

00730 

00731 

00732 

00733 

00734 
00736 

00736 

00737 
00736 

00739 

00740 


45 
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i 


w 


. i 



40 KMU a BESL2I VALUE I 00741 

60 TO 60 00742 

50 KMU • BESL1I VALUE I 00743 

C 00744 

60 IH KMU I 00* 6S* 70 00745 

65 RIGHT a VALUE 00746 

60 TO 130 00747 

70 IPLUS • 1 0074B 

60 TO 90 00749 

80 IM1NUS a 1 00750 

C 00751 

C DtTtRMINt IF LEFT ANO RIGHT BRACKETS HAVE BEEN FOUNO. 00752 

C 00753 

90 IF I IPLUS .EO. 1 .AND. 1M1NUS .EQ. 1 ) GO TO 100 00754 

C 00755 

C BRACKETS NOT FOUNO. RECVLE. 00756 

C 00757 

VALUSV a VALUE 00758 

VALUE « OELTA ♦ VALUE 00759 

60 TO 35 00760 

C 00761 

C BRACKETS FOUND* CALCULATE KMU 0076? 

C 00763 

100 LEFT a VALUSV , 00764 

RIGHT a VALUE 00765 

1FI SIGMA I 110* 120* 110 00766 

110 KMU a FAL2IP1 BESL1* LEFT, RIGHT* .001* 0.0* 50* VY I 00767 

GO TO 130 00768 

120 KMU a FALZIPI BESL2* LEFT, RIGHT, .001, 0.0, 50, YV 1 00769 

130 RETURN 00770 

ENU > 00771 

SUuRUUTINE KOCAL 00772 

C * 00773 

C THIS SUBROUTINE CALCULATES THE CHARACTERISTIC NUMBERS KMU AND QMU 00774 

C 00775 

COMMON /MODES/ MODE 1501, MUI50), IWAVEI50) 00776 

COMMON /KUHU/ KMU( 50) , QMUI50) 00777 

COMMON /BE SSL/ 1S1GN, JSIGN, DELKMU, TOL, N, PI 00778 

COMMON /CNSTNT/ DUM142), NMODES, SIGMA, 0UM2I5) 00779 

REAL KMU, KMUPRM 00780 

C 00761 

DO 130 1«1,NM00ES 00742 

C 00763 

C CALCULATE ORDER FOR BESSEL FUNCTION EVALUATION 00764 

C 00785 

M - IABS( MODE 1 1 ) ) 00786 

IF I M .NE. 0 ) GO TO 10 00787 

ISIGN a 1 00744 

JSlfeN a -1 00769 

6RAKTL a .1 00790 

60 TO 20 00791 

10 ISIGN a MODE II) / M 00792 

JSIGN a ISIGN 00793 
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HKAI1 « HHITNEV AIRCRAFT DIVISION 
St .PANL1B.L4 

bKAK a • N 

IH ISIGN .GE. 0 I 


bRAML • N 0079 A 

IH ISIGN .Gfe. 0 I 60 TO 20 > . I 00795 

i. 00796 

L NfeGAT IVfe ORDER. IF EVEN* SIGN OF BESSEL FUNCTION MILL BE *U IF 000* 00797 

C SIGN OF BESSEL FUNCTION MILL BE -1. 00798 

C 00799 

IF! M/2 I *2 .EQ. N I I SIGN • I OOdOO 

C 00601 

20 NUMMUS ■ MU(I) • 1 00807 

C 00903 

C CALCULATE CHARACTERISTIC NUMBER KMU CORRESPONDING TO NOOEI1I AND HUlIlOObwA 
C THE VALUE OF KMU HILL BE THE MUIlldl ROOT OF THE EQUATION DEFINING 00605 
C THE SYSTEM OF SIMULTANEOUS EQUATIONS 00806 

C 00807 

KMUPRM « 0.0 00608 

00 90 J&l* NUMMUS 00909 

IF I M .EG. 0 .AND. J .EQ. I ) GO TO AO 00810 

CALL KMUCALl BRAKTL* OELKMUt KMUPRM* BRAKTR ) 00611 

bRAKTL > BRAKTR 00812 

AO CONTINUE 00813 

KMUI1I - KMUPRM 0081 A 

C 00915 

C CALCULATE CHARACTERISTIC NUMBER QMU CORRESPONDING TO H0DEI1I AND NUI1I00816 


IF THE HUB/1 IP RATIO IS ZERO* SET QMU TO ZERO AND CONTINUE 


IF I SIGMA 1 
60 GMUI1! 

80 1FI KMU 1 1 1 ) 


60* 60* BO 

GO TO 100 
90* 60* 90 


90 CALL BESJI KMUPRM* M-JSI6N, EMM1* TOL* IER I 
EMM1 > IS1GN * JS1GN • EMM1 

CALL BESJl KMUPRM* M* EMJ* TOL* 1ER2 I 

CALL BESYI KMUPRM* M* EMY* IER3 ) 

EMJ s ISIGN • EMJ 

EMY * ISIGN • EMY K 

CALL BESYI KMUPRM* M-JS1GN* EVM1. 1ERA I 
EVM1 b ISIGN • JS1GN * EVM1 

A « EMM 1 - 1 M • JS’GN • EMJ ) ✓ KMUPRM 

b « EVrtl - 1 M * JS1GN • EMY I / KMUPRM 

SIGMAK b SIGMA * KMUPRM 

CALL bESJI SIGMAK* H'JSIGN* EMM1, TOL* IER5 > 

CALL BESJI SIGMAK* M» EMJ* TUL, IER6 I 

CALL BESYI SIGMAK, M* EMY, 1ER7 > 

CALL BESYI SIGMAK* M-JS1CN* EYM1, IER8 I 
EMM1 b ISIGN * JSIGN * EMM l 

EYM1 b ISIGN * JSIGN • EYM1 

EMJ > ISIGN * EMJ 

EMY b ISIGN • EMY 

C b LMM1 -IN* JSIGN ♦ EMJ I / SIGMAK 

0 > EYM1 * I H • JSIGN • EMY I / SIGMAK 

OMUIII b - | A ♦ C * SIGMA I / I B ♦ 0 • SIGMA I 

100 CONTINUE 
9999 RETURN 


0081T 

00B1B 

00619 

00620 
00921 
00822 
00823 
00626 

00625 

00626 
00927 
00828 
00829 
00630 
00831 
00632 
00933 
00836 

00635 

00636 
00837 
00638 
00839 
00860 
00861 
00862 
0086 3 
00666 
00865 
00666 
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MATT S WHITNEY AIRCRAFT 01 VISION 
SC.PANLIB.L4 


VCR 07/25/77 

9*0 12.S0.00 


ENO 00947 

SUBROUTINE PRINT 00B4B 

C 00*49 

C THIS SUBROUTINE PRINTS INPUT AND CALCULATED VALUES 00*50 

C 00851 

COMMON /CNSTNT/ NMEAS* NPREO, NMODES, SIGMA, B, NX, FREQ, A, 00652 

1 ONEGA 00653 

COMMON /REFS/ XRfcF* RREF* THREE 00B54 

COMMON /MEASUR/ XMI50), RMI50), THETAMI50), BETA I 50), PSIISOI 00655 • 

COMMON /PRkDCT/ XP(50), RPI50), THETAPISO) 00*56 

CCMMON /MOOES/ NODE (50), MUI50I, 1WAVEI50) 00*57 

COMMON /EMUS/ EMU<50,50), EMUPI50,50), EMUPRMI50,50I, IEMPRT 00658 

COMMON /ANCLES/ DEbRAD, RADDEG 00659 

COMMON /OUTPUT/ AMPRI50), PHASERI50I 00B60 

COMMON /REFCON/ REFERS 00*61 

COMMON /NAVENU/ RXI50) 00862 

COMMON /KtiMU/ RMUI50), 0MUI50) 00863 

COMMON /APHlMU/ AMUI50), PH1MUI50), ICHECK 00864 

COMMON /uVlATE/ SloMAX, SI OMAR, SIGMAT, SIGMA*, SlGNAP, SIGAMI 50 1,00665 

I S1GTMI 50), SIGAHC 150), 10EV 00866 

COMMON /OERSUM/ ARNSUMI50), PRNSUMI50) , AXNSUM(SO), PXNSUNI50), 00667 

1 ATNSUMI50), PTNSUMI50) , A6NSUMI50) , PBNSUMI50), 00868 

2 APNSUMI50) , PPNSUMI50) 00869 

COMMON /OERIVS/ DAM0KNI5G,50), OAMRNSI 50,50), DPHORNI 50,50), , 00670 

1 DPhRNSI50,50), 0AMDXN(50,5O), DAMXNSI 50,50), 00971 

2 DPhOXN (50, 50 ) » DPHXNS(50,50), 0AMDTn(50,50I, 00872 

3 DAMTNSISO, 501 , DPHOTN(50,50), OPHTNSI5O,50l, 00873 

4 0AM08NI50, 50 1, DAMONS! 50, 50), OPHOBNI 50,50), 00874 

5 DPHBNS(50,50t, DAMDPNI 50, 50 I , 0AMPNS(5O,50), 00675 

6 0PH0PN(50,50t , DPHPNSI 50,50 J 00876 

COMMON /MCONP/ XACDMPI50), RAC0MPI50I, TAC0MPI50) , BAC0MPI50), 00977 

1 PAC0MPI50) , XPC0NPI50), RPC0MP450) , TPC0MPI50), 00678 

2 BPC0MPI50), PPC0MPI50) 00879 

DIMENSION AMUOBI50), DEVDBI50) 00680 

COMPLEX RX 00661 

REAL RMU, MX 00882 

C 00993 

C CONVERT INTERNAL UNITS TO OUTPUT UNITS 00884 

C 00885 

THREE • RADDEG * THREE 00686 

C 00867 

DO 20 l»l, NMEAS 00*88 

THETAMI1) « RADDEG * THETAMII) 00669 

PS1I1) » RADDEG • PS1II) 00890 

20 CONTINUE 00891 

C 00*92 

00 25 I«I, NPREO 00893 

THETAP(I) - RADDEG * THETAPtI) 00894 

PHASERII) - RADDEG * PHASER II) 00895 

ANPR(I) « 20. • AL0G10I ANPR(I) / REFERS ) 00996 

25 CONTINUE 00*97 

C 00698 

DO 3C I«i,NMOOES 00699 
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St.PANLIB.L4 9*0 

KXll I » RAO DEC • KXIII 00900 

30 CONTINUE 00901 

C 00902 

DO 35 I « l, NMODES 00903 

PHIMUlI) * RAOOEC • PHIMUll) 00904 

AHU0BI1) ■ 20. • ALOGIOI AMU II) / REFPRS I 00905 

35 CONTINUE 00906 

C 00907 

IF! iDfcV *LE. 0 I CO TO 50 00908 

00 40 lsltNMODES 00909 

OEVOblli * 20. • ALOGIOI 1» ♦ SI GAM (II / AMU II I I 00910 

40 CONTINUE 00911 

C 00912 

DO 45 Ul.NMOOES 00913 

01VS0R • 1. / AMUI1I 00914 

XALUMP 111 * OIVSOR • SORT! XAC0MPI1) ) 0091s 

RACUMPII1 * OIVSOR * SORT I RACOMPII ) ) 00916 

TACUHPI1I « D1VS0R • SURTI TACOHPII) I 00917 

BACOHP 111 » OIVSOR • SORT! BACONPIII I 00918 

PACOMP 111 > OIVSOR • SORT! PACOMP 111 ) 00919 

XPCOMPIl! b SURTI XPCOMPI I ) I 00920 

RPCllMP III B SORT! RPCOMPI 1 1 I 00921 

TPCuMPI 1 1 b SURTI TPCOMPI 1 1 I 00922 

BPuUMF II) b SCR T | DPCOMPI 1 1 I . 00923 

PPtuWlI! • SUhTI PPCOMPIII 1 00924 

45 tONIlNUE ' 00925 

C 00926 

50 DO 55 l=l,NMODtS 00927 

OIVSOR b 1.0 / AMUIl) 00928 

AkUSCMIll b OIVSOR • SURTI ARNSUMIII ) 00929 

AXUSUMtll «■ OIVSOR • SURTI AXMSUMIII I 00930 

ATNSUMIll b OIVSOR • SURTI A1NSUMI1 1 I , 00931 

AtNSUHIll b OIVSOR * SORT! ABNSUMII I I ’ 00932 

APNSUM(l) ■ OIVSOR • SORT | APNSUMIII ) 00933 

PRNSUMI 1 1 * SURTI PRNSUM1II I 00934 

PXllSONI 1 1 b SORT! PXNSUMI1) I 00935 

PTNSOMI 1 1 b SURTI PTNSUMI I I I 00936 

PbNSUMI 1 1 > SORT! PBNSUMI 1 I ) 00937 

PPNSUNI 1 1 « SORTI PPNSUMI 1 1 I 00938 

55 CONllNUt 00939 

00940 

CONVERT ANY NEGATIVE ANGLES TO POSITIVE ANGLES FOR PRINTING 00941 

00942 

CALL ANGPOSI PHIMU. NMODES I 00943 

CALL ANGPOSI PHASER, NPMED I 00944 

00945 

PRINT INPUT VARIABLES 00946 

00947 

WR1TEI6, 90001 00949 

9000 FORMAT I 1H1, T45, •♦♦♦ MODAL CALCULATION COMPUTER PROGRAM 100949 

WRIT El 6, 9001 1 NMEAS, NPREO* NMODES 00950 


9001 FORMATI //, T56, «... INPUT VARIABLES ... *, //, T5* * NUMBER OF MEA00951 
1 SURE ME NT LOCATIONS b i, 12, T5I, ‘NUMBER OF PREDICTION LOCATIONS >00952 





n n n 





PRATT t WHITNEY AIRCRAFT DIVISION VKR 07/25/77 

SC.PANL1B.L4 *.0 12.50.00 

2 S 12. T96. 'NUMBER OF IMOOE.MU) SETS ■ 12 I 00953 

WR1TEI6.9002) 00954 

9002 FORMAT ( //, IX. «... INPUT MOOES ...». //. T5» ■MODE** T14, 00955 

1 'CIRCUMFERENTIAL*. T34. 'RADIAL*. T47, *WAVE • » /. Tlfc. 'MOOE NUMBE00956 
2R'» 134. 'ORDER*. T45* 'INDICATOR* I 00957 

DO 60 1»1,NM0DES 009SB y 

WR1TEI6.9003) I, MODE II). MUII), 1WAVEI1) 00959 i 

9003 FORMA) I 5X. 12* 1IX» 14. 13X. 12. III. 12 I 00960 

60 CONTINUE 00961 

C 00962 

C PRINT REFERENCE VALUES 00963 

C 00964 

MR1TEI6.9004) REFPRS 00965 

9004 FORMAT) ///, IX* *... REFERENCE VALUES ...** //. TS. • REFERENCE PR00966 

1ESSURE * •» E9.4 ) 00967 

009C8 

PRINT TtST GEOMETRY AND CONDITIONS 00969 

. 00970 

WRI1E 16. 9005) SIGMA. B. MX. FREQ« A* OMEGA 00971 

9005 FORMAT I ///• IX. *... TEST f EOMETRV AND CONDITIONS ...** //. 15. 009)2 

l'HUB / TIP RATIO ■ •• F8.3» T42, 'OUTER RAOIUS OF DUCT * *. F8.2» 00973 

2T84, 'AXIAL MACH NUMBER * *. F8.2, /» TS* 'FREQUENCY * *. FB.2, 00974 

3742. 'SPEED OF SOUND « *. F8.2. T84. *f OIAN FREQUENCY » ', 00975 

4F10.2 ) . 00976 

C 00977 

C PRINT CALCULATED NODAL AMPLITUDES ANO PHASES 00978 

C 00979 

WRITEI6.9000) 00980 

WRITE |6 .9006) 00991 

9006 FORMAT) //, T45* '... MODAL AMPLITUDE ANO PHASE CALCULATION ...*» 00982 

1//. IX. '... CALCULATED MOOAL AMPLITUDES ANO PHASES ...** //. T5. 00983 
2 'MODE • * T12* 'CIRCUMFERENTIAL* , T30. *RAOIAL*. T41, 'WAVE*. T47. 00984 

32I6X. 'AMPLITUDE'), T84, 'PHASE*. T98, 'AXIAL WAVE NUMBER*. ' 00985 

4T126, *KMU*. /, T14, 'MODE NUMBER*. T30. *0R0ER*. 739, 'INDICATOR *0098 6 
5, 753, '(PRESSURE)'. T71, *)DB)*. T82. '(DEGREES) *. T98» 'REAL*. 00997 
6T109, 'IMAGINARY', / I 00988 

DD 80 1=1,NH00ES 00989 

WRITE(6,9007) I, MOOE(I), MUII), IWAVEIU, AMU(I), ANUOBII). 00990 

1 PHIMUI1), KXII), RMUII) 00991 

9007 FORMAT I 5X. 12. 4X» 14, UX, 12, 8X, 12, SX , E12.6, 3X, E12.6. 00992 

14I3X.F10.4) ) 00993 

80 CONTINUE 00994 

C 00995 

C PRINT REFERENCE LOCATION VALUES 00996 

C 00997 

WRITE 16,9008 ) XREF, RREF. ThREF 00998 

9008 FORMAT) ///, IX, '... REFERENCE LOCATION ...*• //. T10. *X', T27, 00999 

1*R>, T42, *THETA*» ✓/, 4X, E12.6, 2I5X.E12.6) ) 01000 

C 01001 

C PRINT MEASUREMENT LOCATION VALUES 01002 

C 01003 

IF) 1CHECK .GT. 0 ) GO TO 110 01004 

WRITE 16,9009) 01005 
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9009 FORMAT I ///, IX, «... MEASUREMENT LOCATIONS ...*t //• T5, 'L0CATI0010Q6 


INS T23, »XS T40» *R*. T55, 'THETA', T70. 'AMPLITUDES T09« 0100T 

2* PHASE ', /, T6, 'NUMBERS T73. »(B)S T89* 'IPSIIS / ) 01009 

00 100 1>1,NMEAS 01009 

WRIT fc 16. 9010) I* XMUI. RMIIT, THETAM(l), BETAll), PSI1I) 01010 

9010 FC«MAT( TXt 12. 3X, 3ISX.F12.6lt 5X, E12.6* 5X, F12.6 T 01011 

100 CONTINUE 01012 

C 01013 

C PRINT PREDICTION LOCATION VALUES 01014 

C 01015 

110 NHlTEtb. 90111 01016 

9011 FIRM AT ( ///. IX. '... PREDICTION LOCATIONS //» 75, 'LOCATIONO 1017 

IS T23. 'XS 140. *R*. T55 « 'THETA*. T70. 'AMPLITUDES T09» 01018 

2 'PHASES /. 16. 'NUMBERS T69, '(RESULTANT)'. 136. '(RESULTANT I*. 01019 
3/ ) 01020 

00 120 1 «l .NPREO 01071 

WRIT E (6.9010) 1. XP(I). RP(I)» THETAP(IS AMPR(I). PHASER) I) 01022 

120 CONTINUE 01023 

C 01024 

C PR. NT SENSITIVITY CALCULATION VALUES IF NOT A CHECK CASE 01025 

C 01026 

1FI ICHECK .GT. 0 ) GO TO 250 01027 

WRITE (6. 9000) 0102B 

WRITE (6.9012 1 , 01029 

9012 FCJKMATI //. 145. '... SENSITIVITY COEFFICIENT CALCULATION ...' ) 01030 

IH IDEV .LE. 0 ) GO TO 190 01031 

C 01032 

C PRINT ERROR SOURCE STANOARO OE VI AT ION VALUES 01033 

C 01034 

WHITE 1 6. VO 17 1 SIGN AX. SlGMAR. SlGMAT . SlGMAB. SlGNAP 01035 

9013 FORMAT! ///, IX. *... ERROR SOURCE STANDARD DEVIATIONS ...'• //• 01036 

IT7, 'SIGMA X*. T24. 'SIGMA RS T39, 'SIGMA THETA*. T58, 'SIGMA B*.01037 
2T74, 'SIGMA PS1'. /, 4X» E12.6. 4(5X,E12.6) ) ' 01038 

C 01039 

C PRINT MUDAL STANDARD DEVIATIONS 010*0 

C 01041 

WRIT E (6* 9014 ) 01042 

9014 FORMAT) ///, IX. «... NORMALIZED STANOARO DEVIATIONS DUE TO ALL ER01043 

1RCW SOURCES ...'» //, T5. 01044 

1 'ROUES T 12. 'CIRCUMFERENTIAL*. T30, 'RA01ALS T41. 'WAVES T52. 01045 

2 'NOR RAH ZED AMPLITUDES 180. 'AMPLITUDE', TI05, 01046 

3'PHASES ✓ , Tl4, 'MODE NURDERS T30, 'ORDER*. T39. 01047 

4'1NI>ICA1GRS T57, 'DEVIATKINS T60, 'DEVIATION*. T103. 'DEVIATION'01098 
5, /. 193, 'ID8IS TI03, '(DEGREES IS / ) 01049 

DO 140 Ul.NHOOES 01050 

WR1T:I6, 90151 I, RODE US MU(1I, 1 WAVE 111, SKAMCII). OEVOBdl. 01091 
1 S1GTMI1) 01052 

9015 FORMAT I 5X. 12, 9X, 14, 11X. 12, 8X. 12, IX, 3(11X,E12.6) ) 01053 

1-90 CONTINUE 01054 

C 01055 

C PRINT MODAL STANDARD DEVIATION COMPONENTS 01056 

C 01057 

WR1TEI6. 9016) 01050 
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90X6 FORMAT I ///, IX. «... NORMALIZED STANDARD DEVIATION COMPONENTS (ER01059 
IRON SOURCE DEVIATION TINES RMS SUN OF NORMALIZED INFLUENCE C0EFFIC01060 


2IENTS) ...» ) 0X061 

MR1TEI6.9017) 0X06? 

90X7 FOKMATI //, T5. •MODE*. T26. •AMPLITUDE DUE TO ERROR IN*. T92. 0X063 

1* PHASE DUL TO ERROR IN*. /, TXS. *X*t T27. 'A*. T3T. 'THETA*. 0X064 

2T51, *0 •» T62. *PSI*. T79, »X». T91. 'R't T10X. 'THETA*. TIX5# 0106a 
3*b», 1126, »PSI», / I 01066 

DO 160 I*1,NM00ES 0X067 

WRIT E (6.90Itt I 1. XACOMPII), RACOMPIIX, TACOMPII), SACOMP(I). 0X068 

1 PACOMPU). XPCOMPIII. RPCOMPU). TPCOMPtl). 01069 

2 epcoMPin, ppcompid oioro 

9010 FORMAT! M. 12. 5I1X.EXX.4I. 4X. M1X.E1X.4) ) 01071 

160 CONTINUE 01072 

C 0X073 

C PRINT INFLUENCE COEFFICIENTS 01074 

C 0X07S 

100 MRITEI6.V019) 01076 

9019 FORMAT! ///, IX. *... RMS SUN OF NORMALIZED INFLUENCE COE FFICIENTSO 1077 

X ...» 1 01070 

NR17i.t6.9017) 0X079 

DO 200 I«1.NMQ0ES 01000 

NR1H 1 6.90X9) 1. AXNSUMI i ) » ARNSUMII). ATNSUMI1). ABNSUMI1) . 01001 

1 APNSUMI1 )» PXNSUMII ) . PRNSUMII). PTNSUMII). 01002 

2 PBNSUNO ). PPNSUM(l) 010N3 

200 CONTINUt 0X084 

C 01005 

c print partial derivatives 01006 

C 0X087 

MRlttlb.9020) 01099 

9020 FORMAT! ///, XX. «... INFLUENCE COEFFICIENTS ! PARTIAL DERIVATIVES 0X009 

1) ...» I 01090 

DO 2 40 I«I.NMEAS ’ 0X09X 

WRITE! 6, 9021 ) 1 0109? 

9021 FORMAT! //. T*5, • INFLUENCE COEFFICIENTS FOR MEASUREMENT LOCATION 0X093 

X », 12 ) 01094 

NR1TE(6.90X7) 01095 

00 220 J«1 »N MODES 01096 

NR1TLI6.9O10) J. OAMOXNI I* J). OANDRN!! »J) . DAMDTNI I.J ). 01097 

X 0AM0BN1 I *J ). DAMOPNII.J). OPHOXNII.J), 01090 

2 OpHORN! I . J ). DPHOTNII.J). DPHO0N !!»<))» 01099 

3 DPHOPNII.J) 01100 

220 CONTINUE 0X101 

240 CONTINUE 01102 

C 01103 

C PRINT CHARACTERISTIC E-fUNCTION VALUES IF REQUESTED 01104 

C 01105 

250 IF! ILMPRT .Lt. 0 ) 60 TO 9999 01106 

WRI 7£!6, 9000) 0110 7 

WR11EI6.9022) 01100 

9022 FORMAT! //. T23. «... CHARACTERISTIC E>FUNCT10N VALUES FOR NODAL A01109 

1MPL1 TuOF AND PHASE CALCULATIONS ...* ) OHIO 

WRITE!*. 9023) I I.I»1.NMOOES > 01111 


noon 
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9023 F0HMATI ///, IX* «... MEASUREMENT LOCATIONS //. IX. •L0CATI001U2 

IN*. 168. *M00tS*. //, AX. LSISX.I2I. / I °! l J 3 

00 260 J-l.NMEAS 01UA 

NR1TH6, 90261 J, I EMUI 1 ,J I.I-l.NMODES I 01115 

9026 FORMAT I 6X» 12. 5X. 15I1X.F7.3I I OI Ht 

260 LUN11N0E 01117 

Wklf L (6,9023) I I»l«l,NN00tS I 01113 

9025 FukMATI ///, IX. *... PREDICTION LOCATIONS ...** //. IX, •LOCATIONOlltV 
I*. T6*), *Mt!Dfc S* , //, 9X, 15I6X.I2), / I 01120 

DO 2 Bo J*1,NPRED 011 £l 

NHlUlo. 90261 J, I EMUPI1.J)*1*I.NN00ES I 0U?2 

200 CONTINUt 01121 

9999 RETURN 01126 

END 01125 

SUBROUTINE SIMtQCI A, C. NA, NO, SNCUL ) 01126 

01127 

THIS SUBROUTINE SOLVES A NA X NA SYSTEM OF SIMULTANEOUS EQUATIONS 01128 

H AVI NO COMPLEX COEFFICIENTS USING GAUSSIAN ELIMINATION METHOD. 01129 

OHIO 

COMPLEX AI50.1I, Cll), SAVE. 2ER0 01131 

DATA ZERO / (0.0,0.01 / 01132 

C 01133 

SNCUL • 0.0 01136 

DO 260 1*1, NA . 01135 

c 01136 

C FIND MAXIMUM ELEMENT IN JTH COLUMN, RONS 1*1 10 NA 01137 

c Oil IB 

JZ • I ♦ 1 01139 

1FI 1 - NA I *0, 100, 20 01160 

20 VALMX > CABS! Afl.lt I 01161 

MZ • I 01162 

DO 60 KZ*JZ,NA 01163 

B * CABS* AlK2.II I 01166 

1 FI VALMX - B I 60. 60, 60 01165 

#0 VALMX « B °!!tt 

MZ » KZ 01167 

60 CONTINUE ' 0H6« 


C INTERCHANGE ROM CONTAINING MAXIMUM WITH ITH RON 

C 

DU 60 1K«1,N8 

SAVE * All.IRI 

A ( I, IK I * A(MZ.IK) 

A (HZ, IK I • SAVE 

80 CONTINUE 
C 

C NORMALIZE ITH RON 
C 

100 IF! REAL I All, II I I 160. 120, 1 

120 1H A1MAGI All, II I ) 160, 160, 1 

C 

C ERROR - COEFFICIENT MATRIX IS SINGULAR 

C 
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PRATT & MHITNEV AIRCRAFT DIVISION 
SC.PANL1B.L4 


07/25/77 
12. SO .00 


ISO SNCUL ■ 1.0 

00 ID W« 

UO DC 220 L2*J2.NB 

AII.L2) - Al I ,L2I / All. II 

IM JI - NO » ISO. 2 AO. 240 

180 00 290 N202.NA 

AIN/ »LZ I - AIN2.LZ | - AIN2.II • AII.LXI 

200 CONTINUE 
220 C0N11NUE 
240 CONTINUE 


SOLVE FOR COEFFICIENTS 


260 DO 280 M/*i,NA 
C IN/ I * 2 

280 CONTINUE 

CINAI * A 

NC > N 

11 * I 

00 320 12*1 .NC 
KK * N 

L2 * N 

CIL2 I * A 

00 300 N*l. II 

cun * c 

KR * K 

300 CONTINUE 

11 - 1 

320 CONTINUE 
9999 RETURN 
END 

subroutine solve 


AINA. MB) 
NA - I 
I 


NA 

NA - 12 
AIL2.NBI 


CCL2I 
KK - 1 


CIRRI aPaILX.KKI 


II ♦ 1 


THIS SUbRUUT INE SETS UP ANO SOLVES THE EQUATION SVSTEN ASSOCIATED 
MEASUREMENT LOCATION PARAMETERS 


COMMON /bESSL/ DUM3ISI. PI 
COMMON /REFS/ XREF. RREF* THREE 

COMMON /MEASUR/ XM ISO I . RM|50», THETAHISOl. BETA 1 501* PSIISOI 

Common /mooes/ mode i so ». dumiiiooi 

CCMMDN /EMUS/ EMUISO.SOI. EMUPI50.S0I, EMUPRMI SO. SO I. 1EMPRT 

COMMON / APH1MU/ AMUI50I, PHIMUI50I. 1CHECK 

COMMON /CNSTNT/ NMEAS. NPRED. NMODES. SIGMA. B. DUM2I Al 

COMMON /MAVENO/ RXI50I 

COMMON /LMATRX/ eohso.su 

COMMON /DV1ATE/ 0UM4|lS5t* IDtV 

COMPLEX KX. EXPNT. EQ1. EQIS0.S1I. ANSWER ISO I 


C SET UP COEFFICIENT MATRIX 
C 

DO AO 1*1, NMEAS 
DX * XMI1I - XREF 

DR * RMI1I - RREF 


01165 

01166 
01167 
0116B 

01169 

01170 

01171 

01172 

01173 

01176 

01175 

01176 

01177 
01170 

01179 

01180 
01181 
01182 
01183 
01186 
011*5 
01186 
01167 
01188 

01189 

01190 

01191 

01192 

01193 

01194 

01195 

01196 

01197 
4ITH0119B 

01199 

01200 
01201 
01202 

01203 

01204 

om j 
01206 
01207 
0120 * 

01209 

01210 
01211 
01212 

01213 

01214 

01215 

01216 
01217 
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PRATT & WhITNEV AIRCRAFT DIVISION 
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OTMfcTA 


• THETAMI I I - THREF 


C CALCULATE CHARACTERISTIC E-FUNCTION VALUES AND DERIVATIVES FOR RPRIME 
C 

10 R PAINE « DR / B 

CALL tHUCALt RPRIME » EMUI1*1I* EMUPRN|I»II V l I 
C 

IS DO 10 J>1*NM00ES 

EXPNT > CMPLXI 0.0* REAL! KXIJI I • OX ♦ MOOEIJI * DTNETA 

EQI I * J) ■ EMUIJ*II » CEXPI EXPNT I • EXPI -OX • 

l A1MAGI KXIJI I I 

EQllitJI •EOII.JI 

20 CONI INUt 

SET UP RIGHT HAND SIDE 

tOll.NMOOESbll - BETAIII • CEXPl CMPLXI 0.0* PSIIII I I 
EG1I I*NM0DES+1I ■ EQI1*NM00ES*II 
40 CONTINUE 
C 

C SOLVE EQUATION SVSTEN 
C 

CALL SINEOCI EO* ANSWER* NHEAS* NN00ESM* SNGULR I 
1EI SNbllLR I 60* 80* 60 

C 

C ERROR - SINGULAR MATRIX. TERMINATE EXECUTION 
C 

60 NM1 • RHODES ♦ I 

WAIT E lb * lOOOl I I E01 1 1 * J I » J*1 *NMI l,l*l.NMEAS I 
10««0 FORMAT I //* SR, 'COEFFICIENT MATRIX IS SINGULAR** I /* IX* 
11GG13.6 ) I 

STOP 

C * 

C CALCULATE AMPLITUDE AND PHASE VALUES 
C 

BO GO 100 1>1*NN0DES 

AMUtl) > CABSI ANSWER! I | I 

PHIMUI 1 1 * ATAN2I AIMAGI ANSNERI II I* REAL I ANSMERHI I I 

100 CONTINUE 
C 

9999 RETURN 
l NO 

SUUROUTINE SENSTVI 0. NOIN I 
C 

C THIS SUBROUTINE CALCULATES THE SENSITIVITY COEFFICIENTS ASSOCIATED 
C WITH THE EOUAT ION SVSTEN 
C 

DIMENSION EMUAVGISO.SOl* 1R0WIS0I, ICOLISO! 

CIMPLEX KX, EUl. TERM* ZERO* SUM, 0IM)IM*ND1MI* OET 
REAL KMU 

COMMON /DLR1VS/ OAHDRNlSO*SOI* OAMRNSI SO,BO I, 0PHDRNI90* SOI, 

1 OP HANS ISO* SO ) * OAMDXNl SO* SOI* DAMXNSISO«SOI« 

2 OPHOxNISO.SO), t>PHXNS|SO*SO), OAMOTNI 50*501, 


01218 

01219 

0122C 

01221 

01222 

01223 

01224 

01225 
101226 

0122T 

01229 

01229 

01230 

01231 

01232 

01233 

01234 

01235 

01236 

01237 
01236 

01239 
017*0 

01241 

01242 

01243 

01244 

01245 
0124b 
01747 

01240 

01249 

01250 

01251 

01252 

01253 

01254 
01755 
01.56 
01257 
01250 

01259 

01260 
01261 
01262 
01263 
01764 

01265 

01266 
01267 
01260 

01269 

01270 



PRATT & WHITNEY AIRCRAFT DIVISION 
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VER 07/25/77 

*•0 12.50.00 


C 

c 

c 


3 DANINS 150*50 I* OFNOTNI 50*501, DPHTNSt 50*501* 01271 

A DAMDBNt 50, 501* OANBN&I 50,501 « OPHOBNI 50,501* 01272 

5 DPHONSI 50,50) * OANOPNI 50, 501 • OANPNSI 50*501* 01273 

6 OPHDPN ( 50, 50 ) * OPHPNSI 50* 50) 01274 

COMMON /DfcRSUM/ ARNSUMI 50) , PRNSUM(50), AXNSUMISO), PXNSUMI50I* 01275 

1 ATNSUMC50), PTNSUMI50), ABNSUMI5Q), PBNSUMI50I* 01276 

2 APNSOMI40), PPNSUMI50) 01277 

COMMON /MCONP/ XACUMPI50), KACOHP 1 50) * TACOMPISO), BACOMP1SOI* 01276 

1 PAC0MPI50), XPCUMPI50), RPCOMPISO), TPC0MPI50I* 0)279 

2 tlPCIlltP 150), PPCCMPI50) 012B0 

COMMON /DVIATE/ SUMAX, SIDMaK, Sll.MAT, S1GMAB, S1GMAP, S1CAMI50)*012B1 

1 SIG1MI50), S1G AMC ( 50 ) * I0EV 012B2 

COMMON /CNSTNT/ NMt AS, NPRCD, NMODtS, SIGMA, B* DUNK 4) 012B3 

COMMON /MEASUR/ 00^21150), BE1AI50), PSH 50) 012B4 

COMMON /APHIMII/ AMIII50), PH1MUI50) 01285 

COMMON /EMUS/ EMUl50,5il|, EMUPISO.SO), EMUPKMI 50,501* 1EMPR7 01286 

COMMON /REELUN/ RtlPRS 0I2B7 

COMMON /ANulES/ Ml. RAO, RADDIG 0128E 

COMMON /AUMu/ XMll|50), QNUI50) 01259 

COMM UN /WAVENU/ KX (50) 01290 

COMMUN /LHA1HX/ 101(50,51) 01291 

COMMON /MttOl^/ MODE 150), MU (50) 01292 

DATA ZtKO / 10.0,0.0) / 01293 

, 01294 

CALCULATl INVERSE OF Mt ASURkMENT LOCATION MATRIX 01295 

01296 

DO 10 J. 1 ,N01M 01297 

DO ID 1-1.NU1H 01290 

Dll.J) * E01I1 ,J) 01299 


10 CONI INUt 01300 

CALL iNVtRl ( Q, ND1M, OtT* IRON* ICOL ) 01301 

C 01302 

C CALCULATE AVERAGE CHARACTERISTIC E-FUNCTION VALUES * 01303 

C 01304 

DO 40 J*1,NH00ES 01305 

DO 20 1* 1 .NMEAS 01306 

EMUA VGI I «J ) * KMUIJ) • EMUPRNIJ',1) / I EMUIJ.I) • B ) 01307 


20 C0N11NUE 01306 

40 CUNI1NUE 01309 

C 013)0 

C CALCULATE DERIVATIVES WITH RESPECT TO R 01311 

C 01312 

DO 100 KM,NMEAS 01313 

SUM * ZERO ' 01314 

DO 60 J«1»NM0DES 01315 

SUM * EMUAVGIK, J) • EQllK.J) • AMUlJ) • 01316 

1 CEXPI CMPLXI 0., PHIMUlJ) ) I ♦ SUM 01317 

60 CONTINUE 01316 

DO 60 L-l.NMODES 01319 

TERM - QIL,K) • SUM • CEXPI CMPLXI 0.* -PHIMUIL) ) ) 01320 

damdknik.l) > - reali term ) 01321 

DAMKNS (K ,L ) - OAMLRNIK.L) •* 2 01322 

DPHUKNIK ,L ) ■ - AIMAGI TERM / AMU (l ) ) • R ADO EG 01323 


5 '> 


1 


A 


,1 
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DPHANSIK.LI • OPHORNIK.L) •• 2 01324 

ttO CONTINUE 01325 

100 CONTINUE 01326 

C 01327 

C CALCULATE DERIVATIVES M1TH RESPECT TO X 01328 

C 01329 

OU 160 K»l t NMEAS 01330 

SUN • ZERO 01331 

DO 120 JbI'NMuDES 01332 

SUN - EOltK.Jl • KXlJl * A Hut'll • CEXPt CHPLXl 0.. 01333 

1 PHI HU I J I 1 1 ♦ SUM 01334 

120 CONTINUE 01335 

OU 140 L»ltNNUOES 01336 

TERN * UlLtKI * SUN * CEXPI CHPLXl O.t -PHlHUltl I I 01337 

OAHUXNlKtLl * AINAGI TERM I 01338 

DAMXNSIKtL) * OAMUXNIKtLI ** 2 01339 

DHtUXNIK'L) *• - REAL I TERN / AHUIL) I * RADOEG 01340 

OPHXNSlK.LI * DPHOXNIKtL) *6 2 01341 

140 CONTINUE 01342 

160 CONTINUE 01343 

C 01344 

C CALCULA1 E DERIVATIVES Ml TH RESPECT TO THETA 01345 

C 01346 

DU 220 K-l t NMEAS , 01347 

SUN > ZERO 01348 

do ieo j e i»nnudes 01349 

SUN = EU1 IK«JI * MODE I J I » AHUM) * 01350 

1 CEXPI CHPLXl 0.. PHIHU(J) I I ♦ SUN 01361 

180 CUNUNUE 01352 

DU 200 L*1 tNNODL S 01353 

TERN « OIL »K I • SUM • CEXPI CMPLXI 0., -PHIMUILl I I 01354 

OAMDTNIK,LI * AINAGI TERN I / RADOEG 01355 

DAMI NS IK * L 1 * DAMDTNIKtU ** 2 ’ 01356 

DPHD TNIK f L 1 * -REAL I TERN / ANUILI 1 01357 

DPH1NSIK.L1 - UPHOTNIK'L) ** 2 01358 

200 U.NT1NUE 01359 

220 CUNUNUE 01360 

C 01361 

C CALCULATE DERIVATIVES MI IH RESPECT TO BN 01362 

C 01363 

DU 260 L*1 tNMODES 01364 

DU 240 K*1 »NMfc AS 01365 

term » OIL tK) • CEXPI CHPLXl O.t PS1IRI - PHIMUILl I I 01366 

DAHDbNIKtL) * REAL! TERN I 01367 

DAHbNSlK.Ll - DAHDBN(KtL) *• 2 01368 

DPHUtNIK.LI - AINAGI TERN / ANUILI I • RADOEG 01369 

DPMbNSlK.Ll « OPHDDNIKtL) ♦♦ 2 013T0 

240 CONTINUE 01371 

260 CONTINUE 013T2 

C 01373 

C CALCULATE DERIVATIVES M1TH RESPECT TO PS1 01374 

C 01375 

DU 300 LBltNMODES 01376 
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4. WHITNEY aircraft division 

%» .PANT m.U 


VSR 07/25/77 

9.0 12. $0.00 


00 280 K*1,NMEAS 


URN 

I 

1>AMUPN0C,L| 
OAHPNSIKtL) 
(jPHOPNfk«LI 
DPHPNSIk »L I 
280 CONTINUE 
300 CONTINUE 


> QIL.KI • BFTAIK) * CEXM CNRUI 0., RSItKl 
PH1NUIL) I I 

• - AlNAGI TERN I / RAOOEG 
« OANOPNIKtLI •• 2 

* REAL! TERN / ANUllI > 

« DPHDPNIk.Cl *♦ 2 


C 

C CALCULAl t SONS OF DERIVATIVES 
C 

00 3 AO J*i.NNODES 
SUMH * 0.0 

00 320 1*1.NHEAS 

SUNN s OANRNSUi J) ♦ SUNN 

320 COMlNUt 

ARNSUMI.il ■ SUNN 
340 CONTINUE 
C 

00 380 J*1,NH0UES 
SUNN * 0.0 

00 360 I*1»NMEAS 

SUNN - DPHRNSI1.JI ♦ SUNN 

360 CONTINUE 

PRNSUMIJI * SUNN 

380 CONTINUE 
C 

00 A 20 J*1 tNMUOES 
SUNN * 0.0 

00 AUO l*l,NMcAS 

SUNN B DANXNSI I. J) ♦ SUNN 

AOO CONTINUE 

AXNSUHTJI • SUNN 
A20 CONTINUE 
C 

00 A 60 J*1*NM00ES 
SUNN * 0.0 

00 AAO 1*1 tNNEAS 

SUNN * OPHRNST I* Jl ♦ SUNN 

AAO CONTINUE 

PXNSUNIJ) * SUNN 

AAO CONTINUE 
C 

00 500 J*i«NNODES 
SUNN * 0.0 

DC A80 I-l.NMEAS 

SUNN B DAMlNSTItJI ♦ SUNN 

480 CONTINUE 

AlNSUHl J) * SUNN 

500 CONTINUE 
C 

OU !>*.0 J*1#NN0DES 


01377 

01378 

01379 

01380 
C 138 1 
01382 
01353 
0138A 

01385 

01386 

01387 
01368 
0135V 
013VO 

01391 

01392 

01393 
0139* 

01395 

01396 

01397 
0139b 

01399 

01400 
01A01 
01A02 

01403 

01404 

01405 

01406 

01407 

01408 

01409 

01410 

01411 

01412 

01413 

01414 

01415 

01416 

01417 

01418 

01419 

01420 

01421 

01422 

01423 

01424 

01425 

01426 

01427 

01428 

01429 


vS 




k 
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VER 
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12.40.00 


SUNN - 0.0 0|430 

00 520 I-ltNNEAS 01431 

SUNN ■ DPHTNSIlt J) • SIMM 01432 

520 C0N11NUE 01433 

PTNSUNIJI • SUNN 01434 

540 CONTINUE 01435 

C 01434 

DO 540 J-ltNNOOES 01437 

SUNN - 0.0 01438 

00 540 I*ltNNEAS 01439 

SUNN * DAHBNSt I» Jl ♦ SUNN 01440 

540 CONTINUE 01441 

AbNSUNIJ) ■ SUNN 01442 

580 CUN11NUE 01443 

C 01444 

00 620 J>1 V NNU0ES 01445 

SUNN * 0.0 01446 

DO 600 I*1,NNEAS 01447 

SUNN * DPHBNSIlt.il ♦ SUNN 0144* 

600 CONTINUE 01449 

PbNSUNIJI » SUNN 01450 

620 CONTINUE 01451 

C 01452 

00 660 J«1 tNNODES , 01453 

SUNN « 0.0 01454 

DO 640 1«1«NNEAS 01455 

SUNN * OANPNStlfJI ♦ SUNN 01456 

640 CONTINUE 01457 

APNSUNIJ I « SUNN 0145* 

660 COMINUE 01459 

C 01460 

00 700 3*1 tNNODES 01461 

SUNN * 0.0 * 01462 

00 680 !■! tNNEAS 01463 

SUNN « DPHPNSf I » Jl ♦ SUNN 01464 

690 CONTINUE 01465 

PPNSUMI J ! « SUNN 01466 

700 CONTINUE 01467 

C 0146* 

C CALCULATE COEFFICIENTS OF DEVIATION FOR EACH NUOE IF REQUESTEO 01469 

C 01470 

IF I 1DEV .£0. 0 1 GO TO 9999 01471 

SIuR - S1GNAR ** 2 01472 

SIGX « SIGNAX " 2 01473 

S1GT « SIGNAT *• 2 01474 

SIGB « SlGNAB •• 2 01475 

SIGP - SIGHAP »» 2 01476 

C 01477 

C CALCULATE CONPONENTS OF NODAL STANOARD DEVIATIONS 01479 

C 01479 

DO 720 !>1,NN0DES 01460 

XACOMPIll • AXNSUNI I 1 • SIGX 014«1 

RACONPII) > ARNSUN(I) • SIGR 01482 



5 ‘> 
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TAC0MPI1) « ATNSUNII ) * SIOT 01463 

BACOMP ( J ) ■ A6NSUMI I) * SIG6 01464 

PACOHP (11 » APNSUMII ) • SIGP 01465 

720 CONTINUE 01466 

DO 740 I«l*NMODES 01487 

XPCOMPIII ■ PXNSUMIII * SIGX 01488 

RPCOMPII) « PRNSUMIII * SIGR 01489 

TPCGMPlI) ■ r TNSUMI II * S1GT 01490 

bPCUMP(I) ■ PBNSUMlI) • SIG8 01491 

PPCOMP 1 1 I « PPNSUMll) * SIGP 01492 

740 CONTINUE 014V3 

CO 760 1»I*NM00ES 01494 

S1GAMI1) - SORT! XACOMPt 1 1 ♦ RACONPII) ♦ TACOHPII I ♦ 01495 

1 8ACUNPI 1 I ♦ PACGMPII) ) 01496 

S1GTMII I * SORT! XPCIMPI1I ♦ RPCOMPII) ♦ TPCOMPll) ♦ 01497 

1 BPCOMPl I I ♦ PPC0MPI1) I 01499 

SIGAMCI1) ■ S1GAMI 1 1 / AMU II) 01499 

760 CONTINUE 01600 

C 01501 

9999 RETURN 01502 

tND 01503 

SUBROUTINE INVERT I A* N t D* L, M I 01504 

C 01505 

C THIS SUBROUTINE INVERTS A COMPLEX MATRIX. THIS PROCEDURE HAS AOAPTJEO 01506 

C PROM THE IBM SCIENTIFIC SUBROUTINE PACKAGE 01507 

C 01508 

DIMENSION l lilt Mill 01609 
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